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Abstract. The following is a report on some calculations of the statistical distribu- 
tion of delay times due to a fixed-time traffic signal on a single lane highway. 

In Sec. 1, a model of a traffic light is proposed leading to a set of dynamical equations 
describing a relation between the times at which cars leave the light in terms of the 
times at which they arrive. In See. 2, some equations are derived for the conditional 
probabilities that a car will leave at any specified time if it enters at some given time. 
For this, it is assumed that the time intervals between incoming cars form a set of 
independent random variables and that one seeks only the equilibrium solutions for 
which the arrival time of any individual car has a constant probability density. 

In Sec. 3, a procedure for obtaining approximate solutions of these equations is 
derived which actually gives exact solutions for the special case in which the cars arrive 
at equally spaced time intervals, discussed in Sec. 4. In Secs. 5 and 6 this procedure is 
also applied to obtain first and second approximations in the special case in which cars 
arrive with the maximum disorder in spacing possible for this model. 

It is found that to a first approximation, it makes very little difference what statisti- 
cal assumptions are made if one wishes to calculate the average delay. 

1. Introduction. As an illustration of how one may apply statistical methods to 
the study of traffic problems, we consider a relatively simple model intended to simulate 
the flow of automobiles through a single traffic light on a single lane highway. 

The motions of individual cars are described graphically in Fig. 1. The trajectory 
of each car is represented by plotting its position x as a function of time t. The position 
x = 0 is chosen to be the position of the signal and we represent the state of the traffic 
light by a dark line for those times during which the light is red and by a thin line when 
it is green. (Actually the red and green intervals will be effective red and green intervals 
defined more precisely for the particular model in terms of their influence on the tra- 
jectories of the individual automobiles). 

We make the following as our initial postulate: 

[. The cars approach the traffic light in an ordered sequence, without passing. At 
sufficiently large distance from the light all cars move with the same velocity v» . 

In view of this assumption, the trajectories of Fig. 1 must approach parallel straight 
lines asymptotically for z —~ + © or x— — o. No trajectories will cross anywhere. 


*Received October 4, 1954. 
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Fig. 1. An illustration of a set of trajectories in which the position z is plotted against ¢ for each car, 
numbered 0, 1, 2, --- . Heavy lines represent red intervals of the light at z = 0. Light lines at r = 0 
denote green intervals of the light. 


As a result of I, a somewhat more convenient graphical representation of the tra- 


jectories is obtained if we choose a “time”, r = t — x/v, . r represents graphically a 
coordinate perpendicular to the asymptotes of Fig. 1. If one plots x vs. 7 as in Fig. 2, 
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Fia. 2. The trajectories of Fig. 1 are represented by a plot of z vs. r = ¢ — x/v. The “arrival time”, 

the intercept at z = 0 of the asymptote of the lower part of each trajectory, is represented by 7; . The 

“leaving time’’, the intercept at z = 0 of the asymptote of the upper part of the trajectory, is represented 

by 7.t,%,°°' 34 +7, t+ T -:- denote the possible values of r/ for delayed cars. Car 3 in Figs. 
1 and 2 is not delayed. 


the asymptotes becomes vertical lines. At x = 0, 7 = #. Thus the traffic light itself is 
represented graphically in Fig. 2 exactly as in Fig. 1. Also the values of 7 at which the 
trajectories cross the light in Fig. 2 correspond exactly to the values of ¢ at x = 0 in 
Fig. 1. For x < 0, 7 represents the time at which the car would arrive at the light if it 
were to proceed toward the light at a constant velocity v.. For x > 0, r represents the 
time at which the car would have left the light had it been moving at the velocity v9 
from the time it was at x = 0 until arriving at the position x. Differences in r along the 
trajectory thus measure the delay due to the traffic light; the total delay is simply the 
difference between 7 for xz = +o and forz = —o. 

The detailed shape of the trajectories near the traffic light is of much less practical 
significance than the difference between the initial and final values of 7, the total delay. 
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Rather than attempt to propose detailed dynamical laws of motion for each car, we 
shall propose instead a model which describes only the final values of 7 in terms of the 
initial values and postulates nothing about the intermediate states. 

The dynamics of this problem is thus completely described in terms of the variables 
7; = initial value of 7 for the jth car, 


‘ = final value of 7 for the jth car. 


be | 
I 


Since there is to be no passing, we can number the cars according to the values of r 
so that 
tT < Tha and TT < Tha 


We shall for convenience describe the values of 7; and r/ as the times of arrival at 
and departure from the light respectively and we shall further postulate that: 


IT. Tia. — 7, > 86 Tai — > 6. 


+1 = 


Far in front of or behind the light, cars remain at least a certain minimum distance 
apart. Since all cars move with the same initial and final velocities, we can say equivalently 
that they arrive or leave at times separated by a certain minimum interval of time 6. 

There are certain qualitative features of the motion through the traffic light which 
should be considered as a basis for future postulates that will yield a detailed description 
of the relation between 7 and 7; . 

Imagine that some car, which we shall number as car 0, approaches the light soon 
after it has turned red and that all preceding cars have already cleared the intersection. 
The car comes to a full stop and proceeds after the light turns green. The value of 75 
is determined by the time at which the light turned green and is essentially independent 
of 7, . We extend this notion somewhat and postulate the following: 

IIIa. If the jth car is delayed in its motion either by the car preceding it or the 
traffic light or both and if it is the first car to leave the intersection after some red period, 
then 7’ is independent of 7; . 

Although the postulate conforms to reality in most respects, it is undoubtedly in 
some error when applied to a car that approaches the light just before it turns green. 
Such a car decelerates, but before it comes to a complete stop, the light might turn 
green and it would proceed through the intersection leaving earlier than if it had arrived 
somewhat earlier and been forced to stop completely. Such events undoubtedly will in 
most cases be relatively rare and have little effect upon the overall picture. 

In a similar manner we postulate: 

IIIb. If the jth car is delayed and is the second car to leave the intersection after 
a red period, then r! = r/_, + 6, independent of 7; . Similarly for the third car to leave, 
fourth cars ete. We also assume that if all cars leaving during some green period have 
been delayed, then exactly n cars leave at discreet times r/ independent of the manner 
in which they arrived. 

This assumption also contains some error in that the spacing of the cars leaving the 
light does in practice depend somewhat on whether or not the cars are stopped completely 
by the light or just slowed down. Also n may depend somewhat on the values of 7; . To 
make any more detailed assumptions than these would however make the model ex- 
tremely complicated and probably add very little to the general picture. 
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To complete the model we have yet to consider the car that is delayed neither by 
the light nor the preceding car. 

IV. If 7; occurs during a green period and if 7; > 7/_, + 4, i.e. if the jth car arrives 
at the light at least 6 later than the j7-1th car has left the light, we postulate that r’ = 7; . 
A car does not decelerate unless it is forced to do so either by a red light or its proximity 
to the car in front of it. 

Postulates I to IV give essentially a complete description of the dynamics of our 
model. 

In order to express these postulates symbolically, we introduce the following ad- 
ditional notation (see Fig. 2): 

T = duration of cycle, 

T* = length of the red interval. 

We choose the zero of time at the beginning of some arbitrary red period so that if r 
mod 7 < 7%, the light is red at time r. 

In accordance with postulate IIIa and IIIb we let ¢,; = value of 7’ for the first car 
to leave the light after time 7*, assuming that it was delayed, ¢, = corresponding value 
of 7’ for the second car etc. In general, a delayed car will leave only at the time ¢t; + kT, 
7 =1,---,n,k = some integer, ¢; — ¢;-, = 6. 

The dynamical equations of the system expressed by III and IV specify in effect 
that if car 7 enters at time 7; , it will leave at the earliest possible time thereafter. The 
easiest way of expressing these conditions by equations is to give r/ in terms of 7; and 
7;_, . If we let k/_, denote the cycle in which the j-1th car leaves the light and k, the 
cycle in which the jth car enters so that 


kial <tja<(hjat YT, &T <1 < & + UT, 


v¥j—1 
then the equations of motion are: 

(i) If r; < 7/_, + 6, the jth car is delayed by car 7 — 1. This implies that car 7 — 1 
has also been delayed for some reason since by II, 7; > 7;-, + 6 and therefore r/_, > 
T;-1 . We can therefore write r/_, in the form 

ti-) = ki_,T + ¢, for some I. 
7; is then given by 


r= BT + bias f 241%, 
or 


r= (kKi_,+)7T 4+ 4, if l =n. 
(ii) For 7; > r}_,; + 6 (jth car not delayed by car 7 — 1) then 
T= 7; if 7; mod 7’ > 7*, 
77 = kT +t, if 7; mod 7 < T*. 


(If the light is green, it is not delayed; if the light is red, it leaves in the first position 
after the red.) 

It follows from the above description, that if one knows the trajectory of one car, 
it is a rather simple matter, for any particular values of the entering times 7; , to con- 
struct the trajectories for subsequent sequences of cars. By iteration, one simply sends 
each car in order through free or puts it in the first available “bin” according to the 
conditions left by the preceding car. 
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To try to catalog all possible motions for the entire collection of cars is, however, a 
task of another magnitude. It is a rather uninteresting task also because the questions 
of primary interest are concerned not with the detailed behavior of the system but 
rather with the average behavior resulting from some typical distribution of incoming 
cars. To try to analyze the problem exactly would only yield a catalog of irrelevant 
information which would probably completely obscure the relevant features of the 
system. For this reason it is necessary that we approach the problem now from a statistical 
point of view. 

As a final note on the dynamics of this model, it is interesting to note that the equa- 
tions of motion are not reversible. Although the r/ are determined by the 7; , the converse 
is not true. Many possible sets of 7; will lead to the same 7 . 

2. Some statistical postulates. We shall be concerned here primarily with the 
problem of determining the distribution of delay times r’ — 7, assuming that the incoming 
cars are distributed in such a way that the times between incoming cars 7; — 7;-, form 
a set of independent random variables. Although much of the study of queues of various 
sorts have dealt with incoming flows of the above type [1, 2], in fact usually with a 
Poisson distribution, and have usually considered only questions of delay times or 
lengths of the queues, these are by no means the only problems of practical interest. 
The assumption of statistical independence of differences in arrival times is a reasonable 
one for a long highway free of other obstacles but it is apparent that the flow leaving 
the traffic light is not of this type. In order to study the very interesting problems that 
would arise when the flow from one obstacle encounters a new obstacle, one must not 
only relax the restriction on the type of incoming flow that one considers but one must 
calculate other things than just delay times or queue lengths; in particular, one must 
calculate the statistical distribution of leaving times. This is obviously a problem of 
much greater difficulty. 

To return to the simpler problem, we concentrate our attention on car 0 as some 
arbitrary car arriving at time 7» in the cycle khT’ < to < (ko + 1)T and we seek to 
calculate the probability that this car will leave at a time 7) in the cycle kjT < 75 < 
(ki + 1)T. We let P.(7é | ro) = probability that car 0 leaves at 7; if it enters at time 75 . 
tT) will be a continuous variable, but except for the special case 75 = 7. , 76 will have 
only discreet values kiT + t, . We also define P,(rj | 7.) = probability that car k leaves 
at time 7{ if it enters at time 7, . 

By assuming that (r, — 7,-,) are independent random variables, we guarantee that 
the problem will be a Markov process of order 1. In so doing, we let f(y)dy = probability 
that two consecutive cars arrive at times differing by an amount between y and y + dy. 
The function f completely defines the distribution of incoming cars. Actually we have 
also assumed here that f is independent of k. To do otherwise would be of questionable 
physical interest and a great mathematical inconvenience. We will not at this time 
specify any particular form of the function f but it is to be considered as known and to 


be consistent with assumption II 
fy)=0 for y<é. 


The general procedure of attack on this problem is first to find a series of integral 
equations expressing P,(r{ | 7.) in terms of P,_, (ri-1 | te-1) and f(r, — 7x-1) using the 
dynamical equations. We then seek to find an equilibrium solution, assuming one exists, 
by imposing the additional condition that P,(r{ | 7,) is not a function of k. In this way 
we obtain a set of integral equations involving only P,(7’ | 7) as the unknown. 
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We begin the derivation of the equations by considering P,(r{ | 7,) for r7 = k{T + ¢, 
and / > 1, i.e. we seek to find the probability that car 1 will be delayed but will leave 
at various 7{ other than the first position after some red light if it enters at r, . According 
to the dynamical equations, this implies that car 0 left at 75 = k{T + t,_, , (k6 = kf). 
Thus 


P,(k{T + ty | 71) _ | dro P,(kiT +- ty-1 to) f (11 — To) 
bin (la) 


for % < kT + £., f > k, 


P\(kiT + t, | 71) = 0 for 1 >khT+¢,, > oe (1b) 


[We assume here that 7, can have any value with equal probability. In general one 
obtains an equation of this type for the joint distribution of r{ and 7, in terms of the 
joint distribution of 7{ and 7, . If 7, is uniformly distributed, so also is 7, since f depends 
only upon (7; — 70). Thus one obtains an equation involving only condition probabilities.] 

The upper limit of the integral in (la) is somewhat arbitrary since f(r, — to) = 0 
for 7T > 7; — 6. 

The case / = 1 and the case r{ = 7; must be considered separately. For 1 = 1, car 1 
can leave in a first position k{T + ¢, , if either of the following applies. 


(i) nm < KIT and =o = (ki —DT +h. 
(ii) MT <1 <kT+T* and 6 <T. 


From this we obtain 


P\(kiT + t, | 7) = | | dro Po[(ki — 1)T + t, | tolf(r1 — 7.) for 71 < KT, (2a) 
for iT <7, < kiT + T*, 
PART + i, | 7,;) = 0 for iT +T* <7. (2c) 


The upper limits of these integrals are again somewhat arbitrary. 
Finally, for the case r{ = 7, , we obtain 


P,(7 | 71) = | drof >> P(rh| to) }f(71 — Tr) (3a) 


if 7; moa T > T*. (3b) 
P,(r, | 73) = 0 if tr, mod T < T*. (3c) 


Equations 1, 2 and 3 give a complete description of the function P,(r7{ | 7,) in terms 
of P.(75 | 7o), assuming that 7, is uniformly distributed (consequently also 7,). 

The next step is to seek an equilibrium distribution. This is done simply by dropping 
the subscripts on the function P so that P, is the same function of its argument as P, 
is of its argument. Equations 1, 2 and 3 then become a set of simultaneous linear integral 
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equations. Rather than rewrite these equations without the subscripts, we shall hereafter 
use them with this minor alteration implied. 

As the equations now stand, we have an infinite set of equations, for determining an 
infinite set of functions (one for each of the discrete values of ri) of a contiiuous variable 
1,—-27 <1<7]. 

This can be reduced to more manageable form by taking advantage of the invariance 
of the equations to displacements of the time coordinate by multiples of T. Although 
this invariance does not in itself necessarily guarantee that the only solutions of the 
equations are themselves invariant to such a transformation, i.e. 

Piri +T ln +7) = P(r | 71), 
such solutions are certainly the only ones that would seem to make sense physically. 

In view of this, one need calculate P(r{ | r,) only for 0 < r{ < T, i.e. for ki = 0, but 
for arbitrary 7, . The values of P(r{ | 7,) for other values of r{ are then obtained by 
using the periodicity. 

There is still one other fact that we will use to advantage. Equations 1, 2, 3 are 
linear homogeneous integral equations. The normalization should eventually be chosen 
so that 


> P(rt | n) = 1 
since these are conditional probabilities. By using this, we can eliminate one of the 
unknowns [for example, P(r, | 7,)] thus reducing the number of integral equations by 
one but by so doing making them inhomogeneous. 
3. Analysis of equations. We begin by considering Eq. 1. To simplify notation 


we let 


yi = 1 — AT, 
oe es 
and use the periodicity to obtain 
P(t; | y:) = [ dyPltiuily—yfy for w<t (4a) 
Jo 
= 0 for >t 6> 6. (4b) 


We can iterate this equation by successive substitutions of the left side into the 
right side for appropriate ¢, until we obtain an expression for P(é, | y,) in terms of P(é, | y). 
By so doing, we obtain the following equation which corresponds to the statement: 
if car 0 leaves in the /th position, car (— 7 + 1) must leave in the corresponding first 


position. 


dyP(i|y-—yf"yY wW<t, (5a) 


0 


=0 “> 2 te (5b) 


P(t; | y:) 
in which 
f°Y =f, 


[ dy, f(y — yf" "(w), (6) 


#0 


I 


fy) 


ll 
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f‘?(y) represents simply the probability density for a zeroth car and a jth car to be 
separated by a time interval y. The above integrals are multiple convolutions of the 


function f(y). 
Equation (2a) is quite similar in form to (la) and by the same reasoning gives 


a oO 


P(t, | y,) = | dy P(t, |y + T — yf" for 4, < 0. ( 


/0 


Equation (2b) is less pleasant. By using the normalization condition we can avoid 
considering Eq. (3) and write in (2b) 
f > P(ro| rm} =1— DS Plrh| mm). 
P k's? tT’ o>k’,T7 
The case 75 = 7, is now excluded from the sum because in (2b) we integrate only over 
those 7, for which 
tT <7 < kit + T*. 
But 7, > kiT means 7, > kiT + T™*, since 7) cannot lie in a red interval. Thus 75 > 79 
and so the above sum includes only discrete values of 7/ . 


re) n 
f SY Pilon} =1— DY YPM +t | rm). (8a) 
t’o<k'sT k’o=k’, l=1 
Furthermore since P(r5 | ro) = O for ri < 7. , all terms on the left side vanish if 
kiT <7, ;i.e. if a car enters after k{7' it cannot leave before k{T’. 
f + P(r} r)} =0 17> kT. (Sb) 
eta’ sT 


In terms of the new notation (2b) and (8) give 


P(t, | y,:) = dy I, — > p Pitti | y—-—yn MH for O<y, < T*. (9) 
Jo k=0 L=1 

According to (8b), the quantity in the brackets of (9) must vanish for y, > y. In 
the final solution, Eq. (8b) and (8a) must be self-consistent and we apparently have the 
freedom of using (8b) or not here, as we wish. If we do use (8b), we take y, as the lower 
limit of the integral in (9). This leads to a less elegant form of the integral equation but 
one which is more suitable, at least for approximate calculations. In either case, substi- 
tution of (5) and (7) will reduce (9) to a simple form. 

Equation (9), without (8b) gives 


a 


P(t, | ym) =1- dy P(t |y — ws DS? We, S<9, < 7. (10) 
i=1 } 


Equation (9) with (8b) gives 


ae 


P(t, ¥y,) = | dy f(y) — | dy P(t, 


“vy “Vi 


Y; — Y) po ly, 9), O<n <T*, (11) 


1 
in which 
YU, mw) = fy) 


f°YU, y) = / dy’ f° yfy — y’). (12) 
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Equation (10) or (11) expresses the fact that P(t, | y,) equals the probability that 
no previous car left in position ¢, . In (10) this is expressed as one minus the probability 
that a car, / cars ahead left in position ¢, , summed over / and integrated over all arrival 
times of the /th car ahead of the reference car. In (10) this is expressed as the probability 
that the reference car is the first to arrive after the light turned red, less the probability 
that if all previous cars arrive before the red light, the /th one left in position ¢, ; summed 
over all / and integrated over all arrival times before the red light. 

Finally we rewrite Eq. (2c) in the form 


P(t, | y:) = 0 for y, > T°. (13) 


To summarize, we see that Eq. (5) gives an expression for P(¢, | y,) in terms of P(t, | y) 
for alll > 1, whereas equations (7), (10) or (11) and (13) together form an equation for 
P(t, | y) in terms of itself. The problem is essentially reduced to solving this latter set 
of equations for P(t, | y). 

The set of equations for P(t, | y,) could be incorporated into a single inhomogeneous 
integral equation with a discontinuous kernel, although this would be of doubtful value. 
It might be of some help, however, to substitute (11) into (7). Equation (11) gives 
P(t, | y) for y > 0 in terms of P(t, | y) for y < 0. Substitution of (11) and (13) into (7) 
would then give an equation for P(t, | y) for y < 0 in terms of itself, namely 


nT? 


P(t, | y:) = | dy f(y. + T — y) [ dy’ f(y’) 


“9 


+ [ dy P(t, — ns + y+ T) (14) 


_ | dy f w+ l—-y) VsPuty’, v»} for yi < 0. 
70 i=1 


The above set of integral equations is rather disagreeable in any form. At first one 
might be strongly inclined to express everything in terms of Laplace or Fourier trans- 
forms. This would certainly lead to very simple expressions for the convolutions in Eq. 
(6) and (12) and even for the sums over / in (10), (11) and (14). One also observes that 
Eq. (7) would be invariant to displacements in y, except for the fact that (7) is valid 
only for y, < 0, a condition that is not invariant to such displacements. Similarly 
Eq. (10) and (13) “almost” have this translational symmetry. These supplementary 
conditions are just enough to make this approach rather unpleasant and as yet no way 
has been found to overcome certain difficulties. The other forms of the equations, namely 
(11) and (14) show this lack of symmetry more clearly. 

It seems quite unlikely that one will be able to obtain even a formal closed form 
soiution of these equations except possibly for special functions f(y). One can, however, 
obtain a useful approximation procedure based upon the assumption that the average 
number of incoming cars per cycle is small compared with the critical number n. 

In the limit zero flux, no interaction between cars, P(t, | y,) must obviously be 
either one or zero accordingly as y is in the interval 0 < y < 7™ or not. In the above 
formulas, this approximation obtains when the integrals in (7) and (10) vanish. A 
much better approximation for low flows is obtained, however, from (7) and (11). 
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We take as a first approximation* 
P(t, | 4) ~ Pil | y) = 0 for m <0, (15a) 
i.e. there is zero probability that a car entering in one cycle will leave in a later cycle. 


Substitution of this into (11) then gives as a first approximation for y, > 0 


ay 


+ 


dy f(y) =1-— | dy fly) for O<y,<T*  (15b) 


ll 


P(t, Yi) 


#0 


= 0 y, > T*. 


This first approximation corresponds to the statement that a car leaves in a position 
t, if it arrives during the corresponding red period and if it was also the first car to arrive 
in this period. 

One can make a sequence of approximations based upon this first approximation by 
successive substitution of each new approximation back into the equations to obtain 
the next approximation. Each new approximation will give a correction due to the 
influence of cars arriving in cycles further and further removed from the reference car. 

If one expresses all the P’s in terms of the solution of (14), it suffices to apply this 
procedure simply to this equation. Thus we obtain for y, < 0 


P,(t, | y.) = 0, 
a T* . 


Pt, | ys) = dy f(y. +T — y) | dy’ f(y’), 


Jo Jy 


ss { (16) 
Pit | m1) = dy Plt | —wif YU tm + 7) 
»T* 0 ) 
- | a’ f-mwHtT—-yX KfPyt+y’, y')f for a> 2, 
J0 l=] 
and 
P(t, | ys) = DO Pilts | y)- (17) 
i=1 


Successive corrections generally become increasingly more difficult to evaluate and 
so the value of this procedure will be limited in practice to situations in which the first 
few approximations already give results of sufficient accuracy. 

4. Ordered flow. To illustrate the above scheme, we consider two extreme cases. 
The simplest case is that of a completely ordered arrangement of incoming cars, uniformly 
spaced in time; to be considered in this section. The second case, to be considered in the 
next section represents the opposite extreme, namely that with the maximum disorder 
possible for this model. 

In the former case, d is chosen to be the time interval between consecutive equally 
spaced cars and f(y) is chosen to be 


fly) = iy — d), 6 = Dirac 6-function (18) 
f(y) = oy — ld). 


*The subscripts on P will be used here to denote successive corrections in this approximation scheme 
and are not to be confused with those used in Sec. 1 (and later discarded). 
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It follows immediately from (16) that if nd > T, i.e. if the incoming flow is less than 
the critical flow, then P,(t, | y,) = 0 for y, < 0 and consequently all higher corrections 
also vanish. 

P(t, | y,) is given exactly from (11) by 

P(t; | y,) = 1 if 0 <y, < dand 7* (19a) 
0 otherwise. 


From (5) one then finds 
= (1 if (l—ld<y<ld and y, <t, 
Plt, | y,) = ¢ 
lo otherwise for 6 < Tf", (19b) 


_ fi if (l—ld<y,<(l—ld+T* and y, <4, 
lo otherwise o> Tt. 


Thus the complete solution is in this case quite elementary. We can pursue the 
analysis further and investigate the distribution of delay times. We define p(é) such 
that p(t)dt is the probability of a delay of ¢ to ¢t plus dt. With the uniform distribution of 
arrival time for the reference car, in accordance with earlier postulates, one finds that 


pt) = 7 Pl |i — 9. (20) 


There will always be a non-zero probability for no delay. p(¢) is therefore defined as 
above only for ¢ > 0. The probability of zero delay is related to P(r | +) which was elimi- 
nated from the integral equations by the normalization condition. The probability of 
zero delay is obtained also from a normalization condition and is given by 


1 — [ dt. 


) 
Substitution of (19) into (20) gives first of all 
p(t) = 0 for t<0 or > &< 


For 0 < ¢ < ¢, , there are numerous special cases depending upon the relative size of 
d, T* and ¢, but in any case p(¢) usually has, for various ¢, values of 0, 1 or 2 depending 
upon whether 0, 1 or 2 of the step functions in (19) overlap for that ¢ when substituted 
into (20). In any particular case it is a simple task to find the appropriate explicit form. 
In a typical case p(t) will oscillate between the values 1 and 2 in this interval as illus- 
trated in Fig. 3 for a particular choice of parameters. 

From p(t) one can calculate averages of any function of ¢. In particular, one can 
evaluate the average delay 


(i) = 7 > | t dt P(t, | t, — 2). (21) 


For any given set of constants for the traffic light itself, one can easily calculate (é) 
as a function of d or as a function of 8 = T/nd, the ratio of the flow rate to the critical 
rate. Figure 4 shows the results of an exact calculation for a particular choice of 7%, 
6 and ¢, . Although the graph appears on this scale to be quite smooth, it really has a 
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Fia. 3. Distribution of delay times p(t) plotted vs. ¢ for a traffic light with 7 = 1, T* = 0.4, 6 = 1/20, 
t; = 0.5, n = 10 and a flow parameter 8 = 2/3 (a = 10, d = 1/20). 


- — —- -- ordered flow 
———- disordered flow, first approximation 
— — — - disordered flow, second approximation. 


large number of discontinuities in the derivatives resulting from the numerous dis- 
continuities in the function p(é). 

As a practical procedure for obtaining reasonable estimates of the average delay 
for equally spaced cars, the above method gives much more detail and requires much 
more work than is justified by the crudeness of the model. In most cases, particularly 
when 7 is not a small integer, the above formulation will give results differing only 
slightly from the formulas for (t) given by Clayton [3] on the basis of somewhat simpler 
assumptions. Clayton calculates (¢) on the basis of a completely specified sequence of 
arrival times whereas it has been assumed here that any given car arrives at a randomly 
distributed time but that the differences in arrival times of consecutive cars are specified. 

In terms of the notation of this paper, Clayton’s formula reads 
= (t, — 6/2)° =P (22) 
2T(1 — n 68/T) par 

The above formula is also plotted in Fig. 4 for comparison, using the same parameters. 
It can best be described perhaps as a smooth out version of the previous results. The 
differences between the two curves would however be more accentuated if we had chosen 
to compare them for a light described by a smaller value of n or one for which 6 was not 





(t) = 


so small compared with ¢, . 
5. Disordered flow. The postulates made in Sec. 1 exclude the possibility of applying 


the above procedure directly to a flow which is completely disordered, i.e. one for which 
the spacing of incoming cars obeys a Poisson distribution. Although this type of flow 
has been investigated by several people, it is not obvious that it is representative of the 
true situation especially for high density of cars near critical flow, the only circumstance 
in which the statistical distribution of cars may have appreciable consequences. 

The assumption has been made in Sec. 1 that cars cannot arrive at time intervals 








1956) FLOW OF HIGHWAY TRAFFIC 365 


less than 6. This assumption has been used at several places in the derivation and, at 
least for certain single lane highways, is perhaps a reasonable one even for high density 
flows. 

Instead of considering a Poisson distribution which corresponds to a distribution of 
arrival times with maximum disorder consistent with a known value of the average flow 
or average time interval between cars, we shall consider here a distribution with maximum 
disorder consistent with the minimum separation assumption and a known value of the 
average flow or time interval between cars. 

The distribution satisfying these specifications must be of the form 


—a(t—3s) 


Sly) = ae for y>é (23) 
= 0 for y<6 
in which @ is to be found from the known value of the average of y, denoted again by d 


as in the previous case 


d= 6+ I/a, a =(d— 6)"; (24) 


f(y) is essentially a Poisson distribution with a displaced origin. 
One easily calculates from (6) that 


-a(y-—id) 
ee ; @ 
foi) = ye senses ° 
f(y) a G- Dt (y — 76) for y>jé (25) 
= for y <jé 
and from (12) 
LY, y) = fy) for y<6,j>1 (26) 


per re —%y + 5) for Y; > 5,j > :. 
The sequence of approximations represented by Eqs. (16) and (17) does not yield 
a simple exact solution as it did for the ordered flow, but it does give a very rapidly 
convergent series of approximations except for flows very close to the critical value. 
From the first approximation in (16), one obtains 
— for b’<y, < T* 


P,(t, | y,) = 1 for 0<y<6 (27) 


lo for y¥, <0 or yi >T* 
and for 1 > 2, P,(t; | y,) can be found in terms of the function* 
\ wipedgate 
I{u, p] = rp + 1) [ ve" dv for u>0O (28) 
= 0 for u<0O. 
To calculate p(t), we again use Eq. (20) but it is convenient to modify the approxima- 
tion scheme slightly. Substitution of (5) into (20) gives 


p(t) = = {Pt l4—-—)+ > | dy P(t; |t: —t— ws}. 
l=2 D 


( 





*This function has been tabulated for various values of u and p in Tables of the incomplete T-function 
edited by K. Pearson, London, 1922. 
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If one substitutes (25) into this, one obtains 


er as we (az) *) 
p(t) = 1 {ou a= o> [ adz P(t, |t, — t — ze rit 
\ /0 l=2 — 4):) 
The summation now represents the first (x — 2) terms in the expansion of e“*. It is 


convenient to write 
e~** ¥ (az)? /(l — 2)! = 1 — Tfaz(n — 1)7'?, n — 2). 


If now one calculates a sequence of approximate values of p(t) from the series of 
approximate values of P(t, | y,), the contribution to p(é) arising from J[az(n — 1)~'”’, 
n — 2] is usually small in each such approximation and can logically be incorporated 
at each step into the next higher approximation. Thus we calculate approximate values 


of p(t) according to the following scheme. 


1 f 
pil) = p Plt t, — i) + | adzP,(t, | t, — t—2)), 
. ) 
l 
pit) = 7 4-| a dz P,_,(t, | t; — t — 2] lae(n — 1)”, n — 2] 
29) 


p(t) = D pill 


Substitution of (27) into (29) for the first approximation gives only elementary 


integrations resulting in the expression 


0 for Bo or t<0 

p(t) = |? (1+ alti, — 0} for éi&—-—d<t<t, (30) 
T {1+ a6} for i —-T* <t<i,-— 6 
T {1 +06 —e*"”?} for O<t<t, — T* 


This function is plotted in Fig. 3 where it is compared with the corresponding curve 


for ordered flow. 
The evaluation of the first approximation to (¢) from p,(é) is also elementary and gives 


. (4— ‘ne T+ ad 6 
_ / sie ° i aa a co) —_ 2 _— 4 2 
t), 1 + 20) op OT ¢ oT t; 3): 1) 
In most practical cases, this is also very close to Clayton’s formula. If we assume in 
(31) that 6 and ¢; — 7* are both small compared with ¢, , (31) gives 


(t), ~ (1 + @d)t{/2T 


and if we again let 
Yi ‘4 


, = nd nla’ + §)’ 
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we obtain 
2 
t 


oe Sn 
th ~ oT — n0B/T) 


This is the same as Clayton’s formula (22), in which the same approximation of neglecting 


5 as compared with ¢, is made. 
Equation (31) is plotted in Fig. 4 as a function of 8 along with the corresponding 


results from the last section. 








{ 1 l 1 l j 
0 2 4 6 2 10 





Fic. 4. The average delay (¢) is plotted as a function of the dimensionless traffic volume 8 for the same 

traffic light constants as in Fig. 3. The broken line is Clayton’s formula. Solid curve 1, is the curve for 

ordered flow (Sec. 4). Curve 2 is the first approximation for the disordered flow (Sec. 5) and Curve 3 
is the second approximation for disordered flow (Sec. 6). 


6. Second approximation—disordered flow. From Eq. (16) and those of the last 
section, one obtains the following expression for the second approximation P,(¢, | y;) 


for yi < 0. 
= Ifa(y, + T — né)n"'?,n — 1] 
— Ifa(y, + T — (n + 1)d)(n + 1)°”’, n] 

| for y, < —T + T* + nd 
P(t, | ys)) = Tla(y, + T — nd)n-”?,n — 1] (32) 
| — Toy + T — m+ 1)8)(n + 1)°?, 0] 

— exp [—a(y: + T — (n + 1)d]a*(y, + T — T* — nd)"/n! 
lL for —-T+T*+ni<y, <0. 
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The substitution of (32) into the various other formulas will lead to some rather 
cumbersome expressions. Fortunately, the term causing most of the formal complica- 
tions, namely the term of (32) with the exponential, can usually be neglected. 

If one discards this small term, (32) becomes 
PAt, | y:) = Tla(y, + T — ni)n™'?,n — 1] (32a) 

— Ila(y, + T — (n+ 16)(n + 1) M2 n] for y, < 0. 


To obtain P,(t; | y;) for y; > 0, one must substitute this into (11) and take only the 
terms not already included in the first approximation (27). Such a calculation gives 


Pit, | y:) = —Iley, + T-—- a+ )dn+ 1° %,n) for 0<m < 4, 
= —e *~"Tla(T — né)(n + 1)7'”, n] for $< 9, < T*, (2d) 
= 0 : <i 


From (32a, b), one can now compute the next approximation to P(t, | y,) for! > 1, 
and p,(t), however we shall pass over this somewhat unpleasant and uninteresting phase 
of the calculation. 

The evaluation of (¢). , the second order term for the average delay can be evaluated 
exactly from (33) with no trouble, but we give below only a much simplified approximate 
expression which is correct to within a fractional error of (a6)*/2n or a(t, — T)* X 
exp [— a(7* — 6)]/2n. In the example considered for illustration, the error is actually 
less than 3% of the second order term alone. This expression is 

na(T—nb 
(t), ~a (T — né)(1 + ad) | dy I[yn"'’*,n — 1] (33) 


which can be evaluated numerically with the aid of the relation 


a 


| dy Ifyn™'?,n — 1] = 2lfen"'?, n — 1] — nI[e(n + 1)°'”, nJ. 


(t), + (t)2 is plotted as Curve 3 in Fig. 4 for comparison with previous results. 

7. Conclusions. Even though the analysis so far does not include an estimate of 
the errors arising from the mathematical approximations (this can be done, however), 
certain qualitative conclusions are already quite apparent. 

Figure 3 shows that delays for the ordered flow, the first approximation to the dis- 
ordered flow, and Clayton’s formula differ by an amount which would in most cases be 
considered as negligible in view of the crudeness of the model. Clayton’s formula, being 
the simplest to evaluate, still remains the most useful expression at least for low density 
flows. 

It is apparent that the approximation scheme must break down for 8 — | and that 
{t) must become infinite. Even the second approximation gives a correction that is quite 


sensitive to changes in 6 near 6 = 1 and is very small over a considerable range of £, 
being of order 6” for small 8. The third approximation will be even more sensitive to 6 
near 8 = 1, will be of order 8” for small 8 and will be negligible over a larger range of 


8 than the second approximation. 
The assumption that the cars approach the light with a certain minimum separation, 
has a considerable influence on the behavior of (¢) near 8 = 1. As cars become more densely 
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packed, keeping a certain minimum separation, the uncertainty in the position of the 
cars decreases at a much faster rate than if cars were allowed to have any separation 
with the usual Poisson distribution. The delays calculated here will lie somewhere between 
those one would calculate for the Poisson distribution and those for a completely ordered 
distribution. 

Which type of flow is most representative of the true situation is a difficult question 
to answer. A completely ordered flow is certainly unrealistic even though it will in many 
problems, such as that considered here, give simple, qualitatively correct estimates of 
certain features of the traffic. For low volumes of traffic, the difference between the two 
possible types of disordered flow is small, but, in any case, as far as delays are concerned, 
it makes little difference what kind of flow pattern exists: 

For high density flows, the degree of disorder does make a difference. Quite likely, 
the disordered flow considered here is more realistic for single-lane highways and the 
completely disordered flow more realistic for a multiple-lane highway although the 
dynamics of the latter situation is not too well defined. A real experimental test of this, 
however, would not be possible except for nearly critical flows and since the delay is very 
sensitive to the density, it would be necessary that one have an accurately defined volume 
of traffic in an equilibrium situation with no transients or disturbances of any kind not 
specifically taken into consideration here. This is particularly important for high volumes 
because any slight disturbance of the equilibrium will cause queues to build up very 
rapidly and dissipate very slowly. 

One result of these calculations does seem rather surprising, namely that even to a 
second approximation, in the example considered, the average delay is less than a third 
of a cycle. It is also apparent that one will not obtain average delays of more than about 
one cycle on the basis of this model until the volume of flow has reached a value ex- 
tremely close to the critical value. 

Acknowledgment. ‘The author is indebted to Professor William Prager for instigating 
this study of traffic problems and for performing himself most of the unpleasant tasks 
necessary to get such an investigation started. 
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Advanced mathematics for engineers. By H. W. Reddick and F. H. Miller. John Wiley & 
Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1955. xiv + 548 pp. 
$6.50. 


This is the third edition of a very adequate text, first published in 1938 and which presents to 
undergraduate engineering students a wide variety of advanced mathematical topics and techniques. 
This edition, prepared by F. H. Miller, is essentially the same as the second edition. Several new sections 
and topics have been added, but the great bulk of material is unchanged. The problem lists have been 
extensively revised; some new problems have been added and, for some reason, almost all of the formal 
drill problems have been trivially altered so that the answers (all of which are given) are different from 


those of the corresponding problems in the second edition. 
PeTer CHIARULLI 


Lectures on partial differential equations. By I. G. Petrovsky. Interscience Publishers, 
Inc., New York and London, 1954. x + 245 pp. $5.75. 


This book is a translation, by A. Shenitzer, of Petrovsky’s concise introduction to the subject of 
partial differential equations. Chapter I deals with Cauchy’s problem, S. Kowalewsky’s existence 
theorem, characteristics, E. Holmgren’s uniqueness proof for Cauchy’s problem, canonical forms for 
second order linear partial differential equations in one unknown function of two independent variables, 
and canonical forms for systems of linear first order partial differential equations in two independent 
variables. Chapter II is divided into two parts: (a) Cauchy’s problem in the domain of non-analytic 
functions and (b) vibrations. Part (a) deals with the “correct posing’ of Cauchy’s problem, Cauchy’s 
problem for the wave equation in one, two, and three space dimensions and for hyperbolic systems 
in two independent variables, Lorentz transformations, mathematical foundations of special relativity. 
Part (b) is concerned with vibration problems, the so-called ‘‘mixed”’ problems for the wave equation, 
and specially Fourier’s method (expansion in terms of particular solutions obtained by the method of 
separation of variables) for the vibrating string equation. Chapter III is devoted to elliptic equations, 
and covers Laplace’s equation, potential theory, solution of Dirichlet’s problem for a circle by Poisson’s 
integral. The uniqueness of the solution of Dirichlet’s problem is proved by an elementary method 
(not involving Green’s theorem) due to I. I. Privaloy (Mat. Sbornik (1) 32, 464-469 (1925)) and the 
existence of the solution by the Poincare-Perron method of sub- and super-harmonic functions. The 
difference equation method for the approximate solution of the Dirichlet problem is also considered. 
Parabolic equations are discussed briefly in Chapter IV. At the end of each of the last three chapters 
there is a brief but informative survey of related known results. To quote from Professor Courant’s 
foreword to this volume: “It will be highly welcome to English speaking students that Petrovsky’s 
masterly lectures on this important subject are now being made accessible through the present trans- 


lation from the Russian original.” 
J. B. Draz 


Fiinfstellige Tafeln der Kreis-und Hyperbelfunktionen sowie der Funktionen e* und e~* mit 
den natiirlichen Zahlen als Argument. By Keiichi Hayashi. Walter de Gruyter & Co., 
Berlin, 1955. iv + 182 pp. $3.00. 

As far as this reviewer can determine, this is a photo-offset reproduction of the 1921 edition of this 


work, the contents of which are well known and need not be detailed here. 
W. PRaGER 


(Continued on p. 392) 
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ON THE FLEXURAL VIBRATIONS OF CIRCULAR AND ELLIPTICAL PLATES* 


BY 
WILLIE RUSSELL CALLAHAN* 
General Electric Advanced Electronics Center at Cornell University 


In this paper we express R. D. Mindlin’s version of plate flexure equations, which 
take transverse shear and rotary inertia into account, in general orthogonal curvilinear 
coordinates and then we specialize these to polar and elliptical coordinates in order to 
find the frequency equations for the normal modes of vibration of the circular and 
elliptical plates respectively. In particular, we wish to discover those of the eight natural 
boundary conditions for which the normal modes of vibration are expressible in terms 
of product functions. 

In rectangular coordinates, the bending moments (M, , M,), twisting moments 


(M,, = — M,,) and shear (Q, , Q,) are given by the equations: 
a ae st p( 2M + 2), (a) 
Q. = wan 2 + v.), Q = eon( 2 + ¥,), 


where y, , ¥, , and w are plate displacements; D, G, and u are the plate modulus, shear 
and Poisson’s ratio respectively; h is the plate thickness, and k* = 2*/12 is a constant 
for any plate. 

In the case of free vibrations, Mindlin has shown that w, ¥, , and y, can be expressed 
in terms of the three functions w, , w. , and w; by the following equations: 


w=U,+ U2, 


—(¢-1) 2 — 3) 22 4 Ws 

¥, = (o, — 1) _* (o, — 1) a ay? (b) 
- _ 1 ow _ 1) Oz _ Is 

Wy _ (a; 1) oy + (a. 1 ay ar’ 


where w, and w, are components of the displacement perpendicular to the middle plane 


of the plate, and w, is the potential function which gives rise to the twist about the 
normal to the plane of the plate; 


o, = &(S"' + Rd)", o, = 6(S"' + R&)"', 


*Received May 19, 1954; revised manuscript received May 30, 1955. 
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9 


h° eels : ’ 
where R = 12 (coefficient of rotary inertia), 
Dp , ue , 
S — (coefficient of transverse shear), 


~ k’Gh 


*h , , 
5, = = where p and p are the plate density and angular frequency respectively; 


1/2) 


é; = 405((R + S) + [(R — S)’ + 465*]'7}, 
and 
6; = $85{(R + S) — [(R — S)’ + 48‘]'}. 


Mindlin showed further that the w,; are governed by the following three separated 


wave equations 


(V? + &)w; = 0, i = 1, 2,3, (ce) 
where 
ems FO 3 2 ARs a Ss ') 
V = ax a ay’ and 63 = 1 ae ° 


The functions w; are also linked through the following boundary conditions: One 
member of each of the following three products must be specified on the boundary: 


vM; , ¥,M:, ; wh); ; 
where 
w=w,+ wu, ¥: = ¥, cos 0+ y,sin 8, y, = v, cos 6 — y, sin 8, 
Q; = Q. cos 6+ Q, sin 8, M, = M, cos’ 6+ M,sin’ 6 + 2M,, sin 6 cos 6, (d) 
M,, = (M, — M,) sin @ cos 6 + M,,(cos’ 6 — sin’ 8), 


6 being the angle between the normal to the boundary and the z-axis. 

The classical Lagrange theory of plates is a good approximation only when the 
wave length is large in comparison with the thickness of the plate, and this restricts 
the theory generally to low frequency vibrations. The present theory permits extensions 
to moderately high frequency modes, essentially because it includes coupling between 
flexual and shear motions. 

Transforming equations (a), (b), and (c) into general orthogonal curvilinear co- 


ordinates we have: 


w=UW,+ U2, 


Ow, Ow, OW; 
v: as (co; vr 1)h, ry: + (a2 = lh, re: + he an ; 
dw, Ws 


\7 OWe 0 
¥, = (0, — Ihe an + (¢, — I)he oy —h, OE , 
‘ ow Ow. ow 
k’G v ‘1 pet 2 ss 3 
il es 3 + gh, a + h, an |, 


Q: 
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M = D((o, — ayn S88 + wh SB + Prato, wet + Peles ows) 
+ (0. — of Tus + uhs ~ nan (A) 
_ a — av pend 22 au (ay y oe phe (2x) } dws 
a ro aE an Nat) + \ae an + OE On/L\ae at) \f a an 
+ Rylz, y) ‘+ R,,(z, y) 2s) 
si-w — ayl janet g OF OY ad tH ay d2)| (22) 
(1 D| ( olen true *e +A 
7 ou") dw, v1) 
(2u aE an ~+ R,,(z, y) 2 9g + Reale, » 2 an 
_ dx dy (ay (2 oy , dy oz)| (32)" 
+ (os (ns at ae an \ae) + \ae an + aE An/L\ae 
+ R,,(z, y) 2 2+ R,,(z, 92 =) 


aE 


: 0 ny 122 4 se, + Sule, 0 )]) 


aus 











(22) |} a 


h? es : + na? —sz- ~ + L(x, 2 ir e 


Di 


: + L., (2, ye — = 0, += 1, 2,3, (B) 


where: 


L(x, y) = Es a: (ni § 4 + hs = = (i 22) 
[et ie +e aaa) 
P,,(z, y) = Ree | (22) + t) is (2 (22) ] 3 (8 nS) 
+ L(Y + os) Jee 2) 
+ fee [ (2 % 2u(38) — (3H) | 5 (5) 
+ u(y +112 (6%) c 
R,,(z, y) = = mi( 3) E 9 (ii 2u) af os = (i 2) | 


. oy oy 2 (13 2) a 3 (43 20) 
+n | SS (eS) -S (iS 
me dn LOn On * 8j dn On i aj/ J’ 


< 
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le ate Tae 
Siz, y) = ny22 | (22) +9(2% at \ aj) + ae L\ae hig 
fl + ARTS (9) - BI - TRC) 
22 | an) + 3\an) | an aj] ~ an a , 
2 
7 = (2) (24) t = 1 whenj = £ andi = 2 whenj = 7. 
Specializing equations ( 
polar coordinates h, = 1, h. = 1/r 
— (er, — 1) 2 = oe: 1 dws 
y, _— (o; 1) or + (o2 1) + r 06 ? 


G2 — 1dw, _ dws 


1) and (B) into polar and elliptical coordinates we have: for 


w=WUW+U2, 





_ 71 — 1 dw, 
ve = r 00 + r 00 or ’ 
ow OW 1 dw, D 
@, = Kail 6 et" > + 12ue], ~ 
= ee ew _ | ww, | uw 25) 
M, = D(c of or” + r or > 00° + (¢2 I) or” ” r or + r’ 30° 
Ow, 1 dw, 
+0 - “|? or 06 + r Al 
Sia Law, _ 1 dw 1 aw, _ 1 2s 
7s NDC »|} ardor | + (: o| ara@ rr 80 
_1f aw, _ Law, 1 Tu, lh 
2 L er’ r or r oo |)’ 
’w 1 dw 1 d’w i 
— -— 3 —> + 60; = 0, = 1,2, 3. \ 
or” r or for. om ‘ 2,3 (E) 
For elliptical coordinates 
2 2 l 2 
h=h=7 z te a ee 5 
C“(cosh” — — cos’ n) C*(cosh 2§ — cos 2n) 
a J 
w=wW,+U, y= h] (os — 1 ~ ) + (o, — 1) st + aus, 
df on 
) , Ws Ow 
ae h] (o = Ow piano Ss Ow, ~ aus |, 
On On On 
‘ as Ow, OW, Ow; 
Q: = k am), ae + C2 ae an | 
i Ps faw' dw} = C a _ | , IW, a 2, | 
M; = Di (o IY a2 + yu 52 (1 u)| sinh 2¢ aE sin 2 an 
aw, aw, Che, re" aus |} 
2 Vy .. op C2 _ 9 
+ ese 5 (1 1) sinh 2 aE sin 2n aE (F) 


+ (2 — yo 0, 
y2,2 . 
fa’ w Ch, sin On = + sinh 2 ous JT 
Cc 


+ (1 — Weg, ~ 2 
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1 — 2 3° 1 9 
M, = —>* DiiX(e — f2 7. * — c*ri(sin ne + sinh 2¢ aur) | 


+t - |e Ts _ op? _ Qn —2 + sinh 2¢ awe) | 





dE On = 
dw; _ F'ws _ oo (sin aus) | 
a” _ - ow! 4 2k(cosh 2 — cos 2n)w, = 0, ¢ = 1, 2, 3, (G) 


where 2k; = 6,C, C is the semi-focal length of the elliptical plate. 
We will now find the frequency equations for the normal modes of vibrations of the 
circular and elliptical plates which satisfy the following eight boundary conditions: 


CIRCULAR PLATE ELLIPTICAL PLATE 

(1) ¥y =~¥% =w=0 (clamped plate) %=y,=w=0 

2) = w=Q,=0 ve =, = Qu= 0 

(3) ¥ =M,=w=0 ¥y; = M, =w=0 

(4) ¥ =M,.=Q, =0 ¥: = M, =Q, = 0 (H) 
(5) M, = M,3= Q, = 0 (free plate) M, = M, = Q; =0 

(6) Mj = M,, =w=0 M, = M;, =w=0 

(7) M, = w= Q, = 0 M, = vy, = Q =/0 

(83) M_ = y%=w=0 M,=y,=w=0 


when r = r, and — = & respectively. 
It should be remarked that the boundary conditions that we are assuming are 
particular and that other values could be assumed for the quantities involved if desired. 
By assuming that the solutions of equations (Z) and (G) are expressible as product 
solutions we obtain the following respective pairs of ordinary differential equations: 


oes Fo) , 1d 4 (2 — SoG) = 0 





, dr 
(1) 
| aw.(0) ee = 
| ae + m*w,(6) = 0 (¢ = 1, 2, 3) 
eeu + (a — 2q; cos 2n)w,(n) = 
, (J) 


d’w,(é) F 

wy oo (a — 2q; cosh 2¢)w,() = 0 = 1, 2,3) 
where g; = k? and a and m’ are separation constants. The first of equations (I) is Bessel’s 
equation. The first of equations (J) is Mathieu’s equation and the second is Mathieu’s 
modified equation. The solutions of (I) and (J) are of the respective forms: 
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sin mé 


i = 1, 2,3, (K) 


Il 


w,(r) _ J nl 5; re . » w,(6) 
cos mé 


I 
S 
3 


ll 
i) 
3 


Ce,,(€, 9:), a CEm(N, Qi), a 
w(t) = : wn) = ; (L) 
Se,,(é, qs); b = se,n(7, qi); b bin 


| 
a 
3 


where a,, and b,, are characteristic numbers of the Mathieu functions. 
Hence the solutions of (E) and (G) are products of the solutions in (K) and (L) 


respectively. 
We shall assume the following solutions for (E) and (G) throughout in the solution 


of our two problems: 





wy" (r, 0) = AS? J,,(8; , 7) cos m8, i=1,2, 
(M) 
wi" (r, 0) = AS? J,,(6; , 7) sin mé, i = 8, 
wi” (é, n) —_ CC snailf, i) C@onai(n, qi), (a _ Om = Qon+ /3 
wie, 2) = CO°Cemailé, 93) D> ART?(g,) cos (2r + 1)n fori = 1, 2, 
r=0 
(N) 
ws” (€, n) = Ca S@on+1(€, 9:)SCon41(0, qi); (b = Dan = Denes); 


ws (E, 0) = CL? Semmeilé, Qi) >, Barri (g:) sin (2r + 1)n ~~ fori = 3. 


r=0 

By substituting the assumed solutions (M) in the equations (D), making use of the 

boundary conditions (H), we obtain, after some reductions, the following results for the 
circular plate: 


PrRoBLEM 1:y¥, = ¥, = w = Oforr = 1, (clamped plated). 


AP Tn(b: 5%.) + AD Tn(be 1%) +O = 0, 


. m i 


An’(o, — I)S2(5: , 70) + Am (G2 — 1)S2(82 , 70) + Am’ —Im(ds 1%) = 0, (1) 


lo 


AY (a, — 1) = Tm(81 5%) + AP (a2 — 1) = Tn(52 5 Po) + ATS (8s To) = 0. 
0 To 


PROBLEM 2:y, = ¥, = Q, = Oforr = r% 


AM (oe, — I)S2(8; , Po) + AL (og — 1I)T4( 82 , M0) + AD = In( 5 , To) = 0, 


AM (eo, — 1) = Jm(5; 5%) + AP (a, — 1) = J ml5q 5%.) + AP S(b3 , 7) = 0, (2) 


» m 


= Ji (6; ? To) + An Tod n( 5. » To) = Fg r J m( 5s » To) = 0. 


fim F1 
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ProsLeM 3:7, = M,, = w = Oforr="n. 
AM Tn(b: 5%.) + AP I n(b2 , 7) +90 = 


AQ (oy — II(4: 70) + An(or — YIA(b2 » 70) + Am? mm Jn(5s 570) = 0. 
(3) 


1'P'(G, — | J A8, » To) — - J m( 5; 7) | 


0 


+ Aq, eae ates 


-+- Pg E Tbs »%) — AY 7p, mds ’ To) +4 =  T6(B, ; r) | =0 


Nl 


ProspiteM 4:7, = M, = Q,, = Oforr=n". 


Ae m 


1 oid 25: » To) + AM ood 1b: » To) a _- m( 5s » To) = 0, 


(3) m 


AD (a, — ITa(dr 70) + AL (2 — YIm(82 , 70) + AR? = In(da » 70) = 0, 


A‘’(e, — 1) ae Ji(6:; »%) — 2 J m( 5; 7) | (4) 


fim gg OO _ | J (6s » To) = 73S u( 5: “ 


Pg 
+ ad “L (53 ’ To) — » Iu( 3s » To) oS mn J m( 5s ’ ra | =,40. 
ProsLem 5:M, = M,, = Q, = Oforr = rp (free plate). 


An oJ 1 (5: »To.) + AD God t( 52 » To) +0= 0, 


A“ (a; _ | J 16; » To) — = J m( 5; vr) | 


+ A (oe. -_ | a Ri 6» » To) ‘ini a J m( 52 7) | 


AS 1 P 
+o Eze 170) — = Ta( Oa y Po) + + Tul 85 7) | = 0, (5) 
0 0 


2 
10'(o; -~ pf aca >To) + m Ji(6: 5%) — “sy J (5; no) | 
0 0 





+ A? (a, — | 740, » To) + “4 J1( 52 » To) al ra J m( 5, 7) | 
0 0 


AOU a (2 J1(8s » To) — o J (53 ’ r) | = 0. . 
0 ° 
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Prosuem 6: VM, M,.,=w=Oforr=n. 


An’ Tm(5s 570.) + An’ Im(d2 7) +0 = 


AM (a; nat 1 | 7 Ji( 5, » To) = a J m( 6, ir) | 
0 


+ AS, = ne J2( 8%) — 7 Tn( Bo 5 rd | 


+ Se- seta, 2 ~ 2 a (85, To) + us In( 53 in) | =0 0, 
(1) ” Le , um (6) 
Ra (a, — 1) J'7(6; ? To) + r, Ji( 6; » To) 7 = J m( Ox » To) 
0 0 
+ A%(o, — aa yt) + In(b2 5 70) — ume Im 5a .n) | 
AM" — 0 J!(83 >) — a Pay ro) | = 6. 
ProsLEM 7:M, = ¥, = Q, = Oforr=n. 
AM a, T2(81 570) + AP aT (82 7) + An = Tu(8s 50) = 0, 
(1) m (2)/ , mM ? 
A* (o, — 1) - ce J m( 01 . 1) a A; (G2 — a ‘al m( 52 » Ver + . page pa f ® » To) = 0 
AM (a, — hai te) + . eo a oe os Tn(8; , To) ] (7) 
0 
PF hg ee mee th) + - Jit.) = ume Tal 8p 7) | 
ProsLeM 8:M, = y, = w= O0forr=”n. 
AM Tnlb: 57.) + A? J,(82 5%) +0 = 
a); m 2) m 5) vese : 
Al (01 — 1) = Im(81 y 70) + Am’ (G2 — 1) Jn(82 70) + An Ta(ds , 70) = 0, 
0 
A™(a, — 1)| J(8; 7) FE SE 7%) — ym Jnl 8s 5 To) (8) 
To 0 


ae Fh (o; = y| 740, 9 T>) ot - J 18, 9 1) — or e m( 52 rn) | 
0 


9 


+ Az a-9[™ Jt 5s in) — Iulia 7) | = 0. 


By eliminating the constants from equations (1)-(8) we obtain the frequency equa- 
tions for each of the problems. 
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By substituting the assumed solution (N) into equations (F), making use of the 
boundary conditions (H) and simplifying we obtain the following results for the elliptical 


plate: 
ProBLEM I: y; = y, = w = 0 foré = & (clamped plate). 
Ar Ceansi(Eo 5 G1) + A,” Ceon+ilGo , G2) + 0 = 0, 
Ap? (a, — ICeinsilGo  G) + AP?(o2 — 1)Ceinsi(Go 5 92) 
+ A; (2r + 1)Sen+i(Go » Gs) = 0 (I) 
Ar? (eo, — 1)(2r + 1)Ceen+ilGo 5 G1) + Ar? (o2 — I)(2r + 1)CeonsilEo 5 92) 
+ A,” SeinsilGo » qs) = 0 
for every value of r. 
ProBLEM II: y; = y, = Q; = Oforé = &. 
As? (a, — 1)Ceinsilo 5 1) + AL? (o2 — 1)Ceensi(Go , G2) + Ar Seonsi(Go , 95) = 0, 
AS? (a, — Ce nsilbo » G1) + AL (2 — 1)Ceonsi(€o 5 G2) (11) 
+ Af? (2r + 1)Seon+i(€o , 93) = 0, 
Ay? aCeinsilfo G1) + Ar o2CeinsilGo » G2) + Ar(2r + 1)Seon+i(o , Qs) = 0. 


for every value of r. 
By eliminating the constants in (I) and (II) we obtain the frequency equations for 


problems I and II respectively. 
ProsieM III:y, = M,;, = w = Oforé = &. 
Af? Ceonsilgo » G1) + As” CeansilEo , G2) + 0 = 0, 
Ap? (a, — 1)CeinsilGo » 1) + Ar (o2 — ICeinsi(Eo , G2) 
+ A,?(2r + 1)Sean+i(Go , 95) = 0, 
Ar?(o, — 1){—2(2r + 1)[CeinsilGo , G1) — G Ceon+iGo , 9)] sin (2r + 1), 
+ 2FCein+i(Go , 91) cos (2r + 1),} 


+ A; (2 — 1){ —2(2r + 1)Cei,+1(Eo ’ 2) , G Ceon+1(&o ’ 2)] sin (2r + 1), 
(IIT) 


+ 2FCeinsi(E , G2) cos (2r + 1),} 
— Af {[Seth.s(f> » gs) + (2r + 1)?Seon+i(& , 9s) — 2@Sein+1(€ , Ga)] sin (2r + 1), 
_— 2(2r a 1)F S€2n+1(€o ’ qs) cos (2r + 1),} - 0, 


where 





a 


oe u = : “f 
nom = C*(cosh* € — cos’) — C*(cosh 2 — cos 2n) 
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But since the above equations must be independent of 7 we cannot solve this problem. 
It is found that the same thing is true for the remaining five problems of the elliptical 


plate. 
We can thus sum up our conclusions as follows: 


Conclusion I: The problem of finding the frequency equations for the normal modes 
of vibration for a circular plate under the boundary conditions (H) can be solved in 
closed form and expressed in terms of Bessel functions by assuming product solutions 
for Mindlin’s equations (E). 


Conclusion II: The problem of finding the frequency equations for the normal modes 
of vibration of the elliptical plate, under the boundary conditions (H) can be solved 
in closed form and expressed in terms of Mathieu functions by assuming product solu- 
tions for Mindlin’s equations (G), ONLY in the cases when the boundary conditions 
(H) are independent of the bending and twisting moments M; and M,, respectively. 
Thus, the normal modes for the important case of the free elliptical edge do not appear 
to be expressible as product functions in elliptical coordinates. 
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STRESSES IN A CIRCULAR CYLINDER HAVING A SPHERICAL 
CAVITY UNDER TENSION* 


BY 
CHIH-BING LING 


Aeronautical Research Laboratory, Taiwan, China 


Introduction. The stresses in a large tension member having a spherical cavity 
have been investigated by several authors, including the present writer.'~* In dealing 
with such a problem the surface of the member is commonly assumed to be infinitely 
distant from the cavity so that the presence of the surface produces no effect on the 
stresses in the neighborhood of the cavity. However, this assumption is no longer valid 
when the surface of the member is at a finite distance from the cavity.-The considera- 
tion of the effect of the surface near the cavity naturally increases the analytic complica- 
tion of the problem. 

In the present paper a solution for an infinite circular cylinder having a spherical 
cavity under axial tension will be given. The cavity is assumed to be symmetrically 
located within the cylinder so that the theory of symmetrical strain for solids of revolu- 
tion’ can be applied. The solution is obtained by constructing a stress function which 
satisfies the boundary conditions on the surface of the cylinder as well as at its ends. 
The boundary conditions on the surface of the cavity are satisfied by adjusting the 
coefficients of superposition involved in the solution. Here the stress function is necessarily 
a biharmonic function.® The method of solution is first described. Then, as an illustra- 
tion, numerical examples are given for two different radii of the cavity. In particular, 
the maximum stresses in the cylinder are calculated and shown graphically. 

Method of solution. Denote as usual the cylindrical and spherical coordinates of a 
point by (r, 6, 2) and (p, ¢, 6) respectively. For convenience, r, z and p will be regarded 
as dimensionless quantities referring to a typical length a. They are connected with 


each other by 


z= pcos¢, r= psing. (1) 


Consider an infinite circular cylinder of radius a having a symmetrically located 
spherical cavity of radius \a as shown in Fig. 1. The axis of the cylinder will be taken 


*Received April 15, 1954; revised manuscript received Dec. 29, 1954. 
IR. V. Southwell and H. J. Gough, Concentration of stress in the neighborhood of a small spherical 
flaw, etc., Phil. Mag. 1, 71-97 (1926). 

2J. N. Goodier, Concentration of stress around spherical and cylindrical inclusions and flaws, Trans. 
ASME 55, 39-44 (1933). 

3H. Neuber, 7'heory of notch stresses, Edwards Brothers, 1946, English translation, pp. 102-104. 

4‘C. B. Ling and K. L. Yang, On symmetrical strain in solids of revolution in spherical coordinates, 
Trans. ASME 73, A367-371 (1951). Note that in Eq. [32] and [34], p. 370, the factor (2n? + 2n — 3 + ») 
should be read as (2n? — 1 + »). 

54. E. H. Love, Mathematical theory of elasticity, fourth edition, Dover Publications, 1944, pp. 
274-277. 

6A different method of solution in such a case based on Boussinesq’s approach is to find two harmonic 


functions. See footnote 4. 
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as the axis of z and the center of the cavity as the origin. In the absence of the cavity, 
a uniform axial tension of 7 per unit area would be given by the stress function 


Ta’ > : 

= ———, {8or’z + (1 — 2o)2"} (2 

Xo 6(1 + a) t \ ie) \ ) 
where a is Poisson’s ratio. The method of satisfying the boundary conditions when the 
cavity is present is to construct two sets of biharmonic functions each of which gives 
no traction on the surface of the cylinder and at the same time gives no stress at z infinity 
or the ends of the cylinder. These functions are combined linearly and added to xo . The 


LID 














Fic. 1. The cylinder. 


boundary conditions on the surface of the cavity are then satisfied by adjusting the 
coefficients of superposition attached to the functions. The sets of functions are derived 
by differentiation from two biharmonic functions each of which has a singularity at the 
origin. 
Two sets of biharmonic functions. Consider a biharmonic function as follows: 
a l+u 7 a ' ; 
X=5 log +a {Wi(k)I (kr) + Wolk)krl,(kr)} sin kz dk, (3) 
2 Jo 
where J, are modified Bessel functions of the first kind of order n and y, are arbitrary 
functions, and 


uM = COs @. (4) 


The first term of this function in fact represents a center of radial tension at the 
origin. The subsequent integral is added so as to annul the traction on the surface of 
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the cylinder. This function gives the following normal and tangential stresses at r = 1 
or the surface of the cylinder. 


_ifaf (2 1a 5) 3 
[o,] ao | 2 % ar + , or + a2" or Xi =, 


ae ee ee ee ee 
=a aga [PWG RR® — 10) 


| 


+ palk)k{(L — 20)Io(k) + kI,(k)}] cos kez dk, 


_ Cm) _ 3 a s 3° 
Iteehs = »(2{ \"! (Zs iE *"g =) oh]... 


= — at [ RIL + ve) {RTH + 201 — 0)1,(8)}] sin ke dk. (6) 


(1 +2) Jo 


These stresses are annulled throughout the surface of the. cylinder provided that 
by Fourier transforms, 


Wilk) {kI(k) — Ti(k)} + plk)k{(L — 20)Io(k) + K(k} 








2 1 3 _ 
- [ tq > 2) = a+ aa} cos kz dz, 

2 f° 3zsin ke dz 

Vil ) + v(k){kI(k) + 2(1 — o)I,(k)} = >), +2) (6) 
In view of the relation’, 
. 1-3-5 ---(2n —1) f°” cos kz dz 

+: Of ee eee eg “ae 7 
\ ) k 3 (1 + petiente ( ) 


where K, are modified Bessel functions of the second kind of order n, we find 


2 T(k)K i(k) + I(k) Kol’) 








vik) = | PT) — {k + 21 — 0) TK)’ (8) 
be) a 2 lol) Kol) + {h? + 201 — 0) LK y()/k 3 ;' 
a ~ i(k) — (+ 20 — o) Tih) -“ = ae 


With these values of y, and y» , it appears that the integral in (3) becomes divergent 
at the lower limit. However, the divergence can be removed by modifying the integral 
as follows: 


xs = Flog * + a? [ [fy Lolhr) + yalWkrI(kr)} sin ke + fe] dk, (9) 
—i 


“- “0 


where 


)  —_— 4 = io 30 a 
f.i(k) = aaa \2 — o)K.(k) + K } (10) 


7G. N. Watson, Theory of Bessel functions, 2nd edition, 1944, p. 185, Cambridge University Press. 
Cf. the integral considered by Poisson and Malmstén. Also see p. 80 for the definition of K, . 
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lower limit. This modification is 


such that the integral becomes convergent at the 
in the cylinder as well as on the 


permissible for it produces no effect on the stresses 
convergence at the upper limit. It can be shown that the function x, thus obtained gives 
no stress at z infinity or the ends of the cylinder. 

It is obvious that successive differentiation with respect to z gives functions with 
the desired properties on the surface of the cylinder as well as at z infinity, but by 
symmetry odd derivatives must be excluded since they are even in z and cannot enter 
into the required solution. The function x, itself may be included since it gives no un- 
balanced force at the origin. The set of functions is therefore 

oO 0 
Xi Pe te —aee eo 


Oz OZ 


Again, consider a biharmonic function as follows: 


1 =apt+a [ fy,(k)Io(kr) + w,(k)krI,(kr)} cos ke dk. (11) 


“0 


The first term of this function represents a concentrated force at the origin in the direc- 


tion of z. This function gives the following normal and tangential stresses at r = 1. 
(1 — 2o)z 32 ae 
oh = Ga aet — age t [PML - 10) 
+ Wi(k)k{ (1 — 20)Io(k) + kI,(k)}] sin kz dk, 
A2Z—o 3 F 
[ah = —Ga mat qa at | POL® 
+ Wilk kI.(k) + 20 - o)1,(k)}] cos kz di 12) 


These stresses are annulled throughout the surface of the cylinder provided that 


Walk) (kIo(k) — 1i(K)} + yalK)k{L — 20)To(k) + kl, (k)} 
_ A Pit=-se  & t..., 
= ae | (+2) a4 Zye72/ Sin he i. 


V3(K)I(kK) + valk) {kI(k) + 201 — o)I,(k)} 


| 
—_ 572 cos kz dz. 13) 


Similarly, we find in terms of y, and y, , 
v,(k) = {k° — 2(1 — o)(1 — 20)}¥.(k)/k, 
Wi(k) = {¥,(k) + 411 — o)¥2(k)}/k. (14) 


In order to render the integral convergent at the lower limit, the function x, is to be 
modified as follows: 


X.=apt+a’ | [{wWe(k)Lo(kr) + Wulk)krI,(kr)} cos kz 


+ fi(k)r? + f(k)\(1 — 4k’2’)] dk, (15) 
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1956] 
where 
(i a {eK ~. i= K(k) 
J2h) = 0 a 1 ’ 
m(1 + o) k (16) 
, . Al — o)(1 — 2e) J, Pr 
f(k) = 1+ OF {Ks +- i Kn}. 


It can be shown that the derivative of the function x. with respect to z gives no stress 
at z infinity. 
It is equally obvious that odd derivatives of x. with respect to z also gives functions 


with the desired properties. The set of functions is therefore 


o 3° a 
az X2» az* X25 az? X2» 


Since the use of any constant multiplier does not affect the desired properties, we 


may write the two sets of functions as follows: 


1 a 
Wo = Xi, On = Baw, 2 1), 
’ 1 nei 
w= —Osp plate, 20). (17) 


In expressing the functions in terms of spherical coordinates, the following relations 
are useful’. 
1 


a mann _(kp)** 
I,(kr) sin kz = >» (-1) (Qn + 1)! 


~ » (kp)?"*' J2n(2n + 1) (kp)* } 
-I,(kr) sin kz = Do (-1 OS bp, s(n), 
salen me ( (2n + 4 4n+ 1 in+ 5 Pon+i(u) (18) 


Pons i(u), 


are H9s — ! 
= 2s log : + L _ a aS AY 
dz ae p 
rs ¥ (28+ 1)! , 
a2"?! -_ ~ (48 + 1) p™" P.,-1(u) P2,4,(u)} 


and in particular 


Op 
ae - P,(u), 


where P,, are Legendre functions of the first kind of order n. Consequently, we have 


, ] 2 = Qn 2n 2n+1 
Wo = _* log £ Et a Do Pa + yop”) p* Pans iu), 
~ u n=0 (19) 


a ‘gf >» (u) : = 2n 2n 2 n+ \ 
=> 3 te D>. (a2, + 200") p* 'Pansi(u) 


«sp n=0 


*For the first relation, cf. Ex. 63, p. 362, T. M. MacRobert, Spherical harmonics, 2nd edition, 1947, 
Methuen. For the last two relations, see Exs. 21-22, p. 332, E. T. Whittaker and G. N. Watson, Modern 
analysis, 1927, Cambridge University Press; or see formula 36, p. 105, E. W. Hobson, Theory of spherical 
and ellipsoidal harmonics, 1931, Cambridge University Press. 
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and 
wh = —a’P,(u) + a® Do ("Bo + Sop") p*Pansilu), 
~ (20) 
Pe P= a’ P2,.-1 LL) — Poss s+ ( ) + @° y (78 + wet \ pert ( ) 
Wo, = 7 4. 1) p”" a 2s 2sP )p 2n+1\H/) 5 
where 
nm (-—))” sf J ¥i(k) ZmllD\ rere dk 
O26 = (On) 128)! Jo \Qn +1 4n4+ 1f- m” 
an (—1)*** cr? f pa(k) 2nWalk) | onsoe+2 
Boe = 7 ees ke dk, 
(2n)'"(2s + 1)! Jo (Qn 4+ 1 4n + 1 (21) 
on (—1)*"" e ae 
he aes! semen : (ky k2"*2**3 de 
as (4n + 5)(2n + 1)!(2s)! i voli dk, 
bs (—1)"** i 2 
"85, = (k)k?"***** dk 
“2s * (4n + 5)(2n + 1)!(Qs + 1)! [ vA . 
and in particular 
"a = — | {hya(k) + fi) ak, 
"Bo = | {ysl + fa(0)} ak. 
The stress function. Now, construct the required stress function x as follows: 
=xX0+ +T + ¢, 10 ,Wo, 5 Bz.) (22) 


s=0 


where A,, and B,, are coefficients of superposition to be determined from the remaining 
boundary conditions on the surface of the cavity; the factor 7 being introduced to render 
the coefficients dimensionless. The function x» has been given in (2). 

It is evident that the stress function thus constructed gives a uniform axial tension 
of 7 per unit area at both ends of the cylinder and at the same time gives no traction 
on the surface of the cylinder. Its expression in spherical coordinates is readily found 
as follows: 
=s 


*p. eT - 1 
oa tenia: Se 92 Fr } = Mina? mye 
30(1 =< {8P,(u) + 2(1 — 50)P;(u)} 9 Ta’ Ay log ae 





TT 3 = bat 2n—1( Kt) 
—28 -_- a > (Fi 
iin p> (2 fe “} 2 1 4n-—3 ra 3) 


2n 


+ Ta yw (Con + D.,p°)p”” Pade, 


where 


Con = Do ("a2Aae + "B2eBa,); 


(24) 


D., — p>» ("-V¥2,Ac. + *"$2,Bo,). 
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rm . . . . . y 
rhis function gives the following normal and tangential stresses 


1 f d, dt <-—- A . 1 ¥( F) i—# 2) 
= —<(2 — - —_—i— v4, — —=—; — + — — 
Op 3 \ le a) On x a ap be dp op x 








a | Op p p 
9 9 9 
7 Tr +5 : TPA (w) +7 : {2 + DGs +2 io 
(2n + 1)(2n + 2)(2n — 2+ 40) p» __ 4n(2n’ + 3n — 0) } 1 
(4n + 3)p° = 4n — 1 ne (25) 
(2Qn(2n — 1)(2n + 1), 4n(4n — 2)(8n — 4no + 1 — o) 
_ = oe! — Ol ———<— . Don-2 
( o (4n — l)p 
 (4n + 2)(4n + 5 (Qn? — ‘i 1 — a) ho 4 
a rm + 3 dD... P.,.(u), 


21/2 2 
— a=W) ) -ae Hey a (1 3\(, 2 1=# 2),1 
Tp — a Zz (I c) dp + p ou V*x + p ou uM dp + p Ou x 


nie yo | [Qn + 2 
= 270 — °)*PXw) + 70 — #7 If ar Aan 





3 ; : n=1 
(2n + 2)(2n — 2 + 4e) 2(2n? — 1 + o) sa. 1 " 
(4n + 3)p° 2n An. mae l Bon- 2( —2n+1 (26) 
{(2n — Qn + VD) a _ (4n — 2)(3n — - 4no + 1 — @) D 
I p 2n (4n- per: 1)p° 2n-2 
(4n + 5)(4n* + 4n — 1 + 20) on Tass 
in +3 - | bn( Me) « 


Note that in spherical coordinates, 


vi=42(p2), 42 {a — 4’) 2h (27) 
p Op Op p Op 
The remaining boundary conditions to be satisfied are that on the surface of the 
cavity the normal and tangential stresses vanish identically. By inserting p = A into 
(25) and (26) and equating separately the coefficient of each Legendre function or its 
derivative to zero, a system of linear equations is obtained. The system of equations 
may be replaced by the following system 


A‘. = we don tz = 6:.. + 2n(n + 1)(4n — 1)Ci, 
3 = ae 


+ 2(4n + 1)(4n* + 4n® — rn? —n +1—-—0°)Di,, (n>0), (28) 


Bi-2 = 3 aa 61, + (4n + 1)C3, + (2n — 1)(2n + 1)(4n + 3)Di,, (n> 1), 
« ~~ 2 





°Cf. footnote * for the formulas of stresses in spherical coordinates. 
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where 6,,, = 1 or 0, according as m = n or m ¥ n, and 
\! 2n + 2 ta (2n + 2)(2n — 2 + 4c) RB 
Am = anes Aon 1 ee os» 
8 Ns (4n + 3)" 
9 
i = = B, 


= = Gt on™ 
(Qn — (2. | 2n—2 
a (2n_ I(2n + IA C,, 
‘ in’ + 2n—4neo+1—o ~ 
4(2n — 1)(3n — Ino + 1 _ a)" * F 
(4n =— 1)( 1n° i 2n —— Ino + l wes a) 2n-2 


(29) 


a _(4n + 5)" D 
2» (4n + 3)(4n? + 2n — 4no + 1-0) °™”° 


A formal solution by successive approximations is as follows: Write, for n > 0, 


Ah = SAL, BL = DBs”, 30) 


p=0 
where 
(0) I | 
an — wes bo aT sa ~ 6 nn? 
3 , oo aa 
31) 
5 
Be.’ = = ae, ee 
i 3(7 — 50) 
and by iteration, 
ro... & ( ss = a p~-1) 
Ain = 2n(n + 1)(4n 1)C2, 
+ 2(4n + 1)(4n* + 4n® — rr? —n +1-0)Di”, (n > 0), (32) 
Bs) = (4n + 1)032°” + (2n — 1)(2n + 1)(4n + 3) D3" 


in which Ci‘”"” and D{°~” are computed from Aj,?""’ and Bj,’~"’ by means of (29) 
and (24). 

Naturally, the foregoing method of solving the linear equations is valid as long as 
the two series in (30) are both convergent. From physical consideration alone, it seems 
likely that there will be convergence if \ is less than unity. An analytic proof of con- 
vergence may be established by means of the method used by Howland’ and Knight.” 
However, for the sake of brevity no further detail will be given here. The solution will 
be illustrated by numerical examples which follow. 

Numerical examples. Numerical examples will be given for the cases \ = 
The coefficients in (21) are first computed and shown in Table 1, in which the integrals 
* with the aid of tables of Bessel 


+ and }. 


are integrated numerically from Gregory’s formula’ 
functions. Here the Poisson’s ratio o is taken as 0.3. Note that °a,, and °8,, are not 


~R, C. J. Howland, Stresses in a plate containing an infinite row of holes, Proc. Roy. Soc. A148, 
471-491, 1935. 

uR. C. Knight, On the stresses in a perforated strip, Quart. J. Math., Oxford, 5, 255-268, 1934. 

2, T. Whittaker and G. Robinson, Calculus of observations, 4th ed., 1948, Blackie, p. 143. 
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TABLE 1 
| =e eT 
Coeff. 2s 2n = 0 2n=2 | 2n=4 2n=6 | 2#=8 2n = 10 
| ae | | | 
0 —0.50722 | 0.058547|—0.0099477 0.0018554 | —0.00036549 
2 0.90856 | —0.32478 | 0.10476 | —0.031758 0.0092718 
‘a2, 4 . —0.73705 | 0.50010 |—0.26061 0.11653 | —0.047229 
6 “ 0.42356 | —0.46333 | 0.35548 | —0.22024 0.11819 
8 - —0.20102 | 0.32370 |—0.34379 | 0.28190 | —0.19336 
10 -- 0.08479 | —0.18884 | 0.26537 | —0.27813 0.23762 
0 —0.42035 | 0.12257 |—0.032782 | 0.0084192 | —0.0021294 
2 — 0.54584 | —0.31634 | 0.14157 | —0.055231 0.019772 
"Bo, 4 — —0.40438 | 0.38274 |—0.25723 | 0.14117 —0.067883 
6 om 0.22351 | —0.31586 | 0.29768 | —0.21899 | 0.13621 
8 ii —0. 10406 0.20626 |—0.26027 0.24757 | —0.19309 
10 = 0.04351 | —0.11520 | 0.18784 | —0.22390 | 0.21446 
0 | —0.35018 | 0.070237) —0.017677| 0.0044470| —0.0011071 0.00027445 
2 0.37928 | —0.25533 0.12212 |—0.049231 | 0.017970 | —0.0061446 
“28 4 | —0.22980 0.29397 | —0.22533 | 0.13319 | —0.067055 | 0.030298 
6 0.10584 —0.21698 0.24384 |—0.19880 0.13226 —0.076390 
8 | —0.04185 0.12579 | —0.19497 | 0.21005 | —0.17863 | 0.12821 
10 0.01509 | —0.06258 0.12819 |—0.17653 0.18655 | —0.16314 
0 | —0.048364! 0.089523) —0.039753) 0.014241 | —0.0046136 | 0.0014109 
2 0.16114 —0.19141 0.13036 |—0.068389 0.030793 —0.012550 
2n5o, 4 —0.10336 0.18829 | —0.18781 0.13694 | —0.082175 0.043140 
6 0.04842 | —0.12918 0.17907 |—0.17402 0.13451 | —0.088482 
g —0.01938 | 0.07185 —0.13274 | 0.16617 | —0.16093 0.12978 
10 0.00705 —0.03486 0.08297 |—0.13012 0.15450 —0. 14990 
| | 
TABLE 2 
" 2n A$, B3,, Cr. Db, 
0 —3.3487 x 107} 3.5515 XK 107! -— —8.456 xX 10-4 
2 8.4785 X 1071 | —6.037 x 1073 1.144 X 107? —2.411 xX 10° 
3 4 —5.563 xX 107? §.318 xX< 10° —8.070 x 10 7.428 X 10-* 
6 1.072 xX 107? —3.931 x 10-5 5.261 X 107-5 —2.897 XK 1077 
8 1.380 x 10-3 2.669 x 10-6 —3.204 K 10-6 1.267 X 10-8 
10 -. @ & 2S ee ee 1.913 X 107-7 —5.920 XK 107!° 
0 —3.3312 X 107 3.0966 * 107! a +1.158 X 1074 
2 7.4303 X 1077 | —2.283 x 10-4 1.360 * 10-3 —8.133 & 10-¢ 
3 4 —2.125 x 10-3 5.645 xX 10-§ | —2.655 x 10-5 6.424 X 10-8 
6 1.145 X10 | —1.172 x 107 4.604 X 1077 —6.475 X 107° 
8 —4.135 x 10-8 2.240 x 10-* | —7.412 x 107° 7.329 X 10-2 
10 ee oe eee 1.163 * 107% —8.861 X 10-" 
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needed in computation for they are the coefficients of the trivial term pP,(u) or z. The 
coefficients Aj, , Bi, , Ci, and Di, are then computed by the method of successive 
approximations for \ = 34 and } respectively. The results are shown in Table 2. In the 
case \ = 3} the computation has been carried out up to the fifth approximation. The 
convergence appears to be more rapid when X is smaller. The coefficients are readily 
converted to A2, , B2, , C2, and D2, by means of (29). 

It is now rather straightforward to compute the stress at any point in the cylinder. 
The most important one, however, is the maximum stress occurring at the surface of the 
cavity across the minimum section where (p, u) = (A, 0). This stress is given by 


1 (1 —n’ a Pt 9” 9  1l—yw ~a 
max og, = 3 K ns o)( moet -I = —" < \o% + . 2 (. ; ae ees = - 
a p Op Op Op Op p On 


* (. @ 4 . = om a) z. a # 2. | 
Op p On/\pdp- p On/” jr. 


(33) 
1 am wa J (2n + 1)Aon + (Qn — 1 + 20)B,, 
oF ee o7 a ——— L 
— 4nd" “C2, — 2(2n° + 8n — 4n0o + 7 — 5a)A"" Dan pPin+i(0), 

where 

o»\ 2 

P3,.:(0) = (—1(?") = 3 (34) 

n 2 

The following values are obtained. 
N 0, 1 1 


T~* (max os) | 2.045, 2.081, 2.359 
The value in the limiting case \ = 0 is the known result of a large tension member having 
a spherical cavity, that is (27-150) /(14-10c). 














2.5 
K 
? —— ee 
| | 
Sr SS = I Se 
| | 
1 | | 
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Fic. 2. Stress concentration factor K versus X. 
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If we define a stress concentration factor K as the ratio of the maximum stress across 

the minimum section of the cylinder to the mean stress across the same section, then 
K = (1 — d*)T™ (max ¢,). (35) 


Consequently, we have 


] 


tol 


rn | 0, 1, 


J 





K 2.045, 1.951, 1.769, 1 
Here the value in the limiting case \ = 1 can be visualized readily from physical con- 
sideration of the cylinder. Figure 2 shows the results graphically. 

The writer wishes to thank Mr. T. C. Lee for his invaluable assistance in preparing 


the manuscript. 
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BOOK REVIEWS 


(Continued from Pp. 870) 


Transactions of the symposium on computing, mechanics, statistics and partial differential 
equations. Vol. II. Symposium on Applied Mathematics Sponsored by The American 
Mathematical Society and Office of Ordnance Research, U.S. Army. 216 pp. $5.00. 


These transactions of the second symposium on applied mathematics are reprinted from the Com- 
munications on Pure and Applied Mathematics Vol. VIII, No. 1 (1955). The eleven papers presented 
cover a wide range of topics, many in relatively new fields of interest or fields that lie in between well 
established ones. One very interesting feature of the group of papers is the surprising amount of overlap 
in the topics despite the tremendous scope of the conference. Indeed many of the papers are primarily 
descriptive surveys and seem to concentrate attention on the amalgamation of diverse topics into new 
fields. 

The following lists the authors and the general topic or condensed title of the paper: P. M. Morse 
(operations research), J. Neyman (inductive inference), H. O. Hartley (analysis of variance), J. E. 
Mayer (statistical mechanics), M. R. Hestenes (iterative computation methods), J. Todd (numerical 
analysis), A. A. Bennett (ordnance computations), C. A. Truesdell (rate theory of elasticity), J. J. 
Stoker (stability of mechanical systems), F. Bureau (divergent integrals), and W. Feller (differential 


operators). 
G. F. NEWELL 


Engineering analysis. By D. W. Ver Planck and B. R. Teare, Jr. John Wiley & Sons, 
Ine., New York, and Chapman & Hall, Ltd., London, 1954. xii + 344 pp. $6.00. 


This book is intended as an aid for students taking a course in “‘engineering analysis” such as the 
authors have been teaching for many years. The authors recognize that, throughout most engineering 
curricula, instructors try to present their students with carefully constructed and clearly worded problems 
whose solution will serve as exercises in acquiring familiarity with new principles and methods. In pro- 
fessional practice the engineer is never confronted with nicely stated problems, but rather with situations 
involving numerous uncertainties (physical, financial, and human) out of which he must construct his 
own problems and find their solution in a form which enables him to choose a definite plan of action. 
The sort of course which the authors teach is aimed at introducing the student to these facts of engineer- 
ing life by posing him a series of practical situations, stated informally and even somewhat vaguely. 
The student must formulate precise problems, apply the basic principles of mechanics, thermodynamics, 
and electricity and magnetism, set up definite mathematical problems, and finally solve them in a usable 
way. 
This book provides many interesting examples of such situations. It provides also many hints on 
how to attack them in an efficient and orderly way, as well as chapters reviewing the essential tools. 
One chapter summarizes basic physical principles (of mechanics, thermodynamics and heat flow, 
electricity and magnetism, and electric circuits). Another chapter reviews some basic mathematics, 
particularly standard elementary methods of solving ordinary differential equations. A final chapter 
gives useful suggestions on the interpretation and checking of mathematical results. The most interesting 
parts of the book will probably be those in which examples are studied in detail of the breaking down 
of complex situations into solvable mathematical problems whose answers enable the engineer to arrive 
at a decision. This is the most difficult part of professional practice, and one wonders how much can be 
learned by formal study and how much can be learned only from personal experience, more or less 
painful. However, there seems no doubt that the authors have made a valiant and commendable effort 
to help the student realize what sort of dilemmas will confront him in practice, and at the same time help 
him to develop confidence that with common sense, good working habits, and knowledge of some basic 
physical laws, he can make satisfactory progress in analyzing them. 
P. S. SyMoNnps 
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TWO DIMENSIONAL SINK FLOW OF A VISCOUS, HEAT-CONDUCTING, 
COMPRESSIBLE FLUID; CYLINDRICAL SHOCK WAVES* 


BY 
T. YAO-TSU WU 


California Institute of Technology 


Introduction. The problem of the steady cylindrical source-type or sink-type flow 
has been of interest to fluid-dynamicists for several reasons. First, it is known that the 
corresponding problem of an inviscid compressible fluid has an exact solution containing 
a limit line of rather special type, namely, the sonic circle [2]. To the exterior of this 
circle the solution has two branches of values, one has its stagnation point at infinity 
(subsonic branch) and the other starts with maximum velocity at infinity (supersonic 
branch). Both of these two branches terminate at the limit line with infinite velocity 
gradient. Therefore the viscous and heat-conductive effects are expected to play an 
important role in continuing the solution further inward. Second, because of its cylindrical 
symmetry, this problem is one of the few nonlinear flows in more than one dimension 
for which there is only one independent variable, the radial distance. Consequently, 
the equations are simple enough to allow a unified discussion of the various effects. 
These are perhaps the reasons why this problem has attracted the attention of several 
authors [3, 4, 5]. 

In the first part of this investigation the Navier-Stokes equations are given for the 
cylindrical sink flow of a viscous, heat-conducting, perfect gas. Then, with some simpli- 
fying assumptions, the qualitative properties of the solutions are discussed in detail for 
the case of large Reynolds number Re. The definition of Re is Re = p,a,r,/u, , where 
r = r, locates the inviscid sonic circle with sonic speed a, , and p, , u, are the fluid density 
and viscosity coefficient at r = r, . These basic properties of the solutions thus compre- 
hended serve as a useful guide in our final calculation of the solution. 

In the second part of this paper the detailed calculation of the solution is carried out 
for the case of large Re. It is found that there is no single expression available for the 
solution uniformly valid in the entire flow region. The calculation is then performed in 
three different regions characterized by the length r, and the parameter Re. For 
r > r, + O(Re~*”) the approximate solution is obtained by using the PLK method**. 
The result fails to be a good approximation for r too close tor, . For |r — r, | < O(Re~*”*), 
a different similarity rule for the variables leads to a system of cylindrical transonic 
equations which governs the flow across the sonic region. These equations can be inte- 
grated analytically for each of the different order terms. The result shows that the 
solutions belonging to the supersonic branch all contain cylindrical shock-type flow 
in this transonic region. In other words, these solutions gradually deviate away from 
the inviscid supersonic branch, reach a minimum near r = r, and then approach asymp- 
totically to the viscous subsonic branch. It is also found that the shock strength is of 
O(Re~’”*) while the shock-thickness is of O(Re~*”*), results which are quite different 





*Received Sept. 10, 1954. This paper is the slightly condensed version of a Hydrodynamics Labora- 
tory report [1] of California Institute of Technology. The work was sponsored by the Office of Naval 
Research under the Contract Number N6onr-24420 (NR 062-059). 

**This terminology was introduced by Prof. H. 8. Tsien in a series of seminars, held in California 
Institute of Technology in 1954, in which the method due to Poincaré-Lighthill-Kuo was discussed. 
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from those of the plane normal shock. Within this region, the thermodynamic variables 
satisfy the isentropic relation up to the order O(Re™'’*) but deviate from it by a quantity 
of O(Re~*”*). The approximate solution for r < r, — O(Re~*”) is subsequently carried 
out. Finally, the entropy variation of the fluid, and the effect due to variation in viscosity 
coefficients, are discussed. 

The corresponding source flow problem was previously solved, using a numerical 
method, by Sakurai [4]; a qualitative investigation on this problem was later elaborated 
on many points by Levey [5], by making use of some conventional methods in nonlinear 
differential equations. In the latter work, the orders of magnitude of many flow quantities 
of interest were estimated. The present investigation on the sink flow is not merely a 
special case of the cylindrical flow other than the source type, but also presents an 
improved method which is more powerful than those used in the previous works (e.g. 
Ref. [3, 4, 5]). The PLK-method applied to the outer region yields a set of reliable 
boundary values for the transonic region and thus enables all flow quantities of interest 
to be calculated quantitatively in all regions. 

The author thanks Prof. H. S. Tsien for suggesting the problem and Profs. M. 8S. 
Plesset and C. R. DePrima for their assistance on many points. 

1. The fundamental equations. Here we are concerned with the two-dimensional 
sink flow of a viscous, compressible, heat-conducting fluid with polar symmetry. The 
only independent variable is the radial distance r from the origin. The radial velocity, 
u, is the only velocity component and is always negative for sink flow. Let p, p, 7, nu, 
nu’, A, R, C, , C, denote respectively the pressure, density, absolute temperature, co- 
efficients of shear and bulk viscosity, heat conductivity, gas constant, specific heats at 
constant volume and pressure. The momentum equation is 


du _ _ dp al» a eee | ats | 
Mee PSL a ts ~ PFS] Tt Ma it) 


dr r dr 


and the energy equation is 
d (w ” d dT du? 2 (3 du? “) | 
a fe ( = — " eae ll asian = a = see ae 9 
pur dr (3 T 7) dr a dr re dr 3 (u ) 2 dr + r f (1.2) 


The continuity equation can be written in the following form if m denotes the sink 
strength, 


2rpur = —™m. (1.3) 


The equation of state is assumed to be that of a perfect gas, 


p = Rol. (1.4) 


Equations (1.1 — 1.4) are a system of nonlinear differential equations for four variables 
u, p, p and T if uw, »’, \ and C, are known functions of 7’. 
To reduce the equations to nondimensional form, the following nondimensional 
quantities are introduced: 
r* =r/r,, n = logr*, w= —u/a,, 6 = T/T, = (a/a,)’, 
(1.5) 


p* = p/p , p* = p/p. , u* = p/m, u* = p/m, 
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where quantities with the subscript 1 are fictitious quantities which would occur at the 
local Mach number unity for nonviscous and non-heat-conducting gas. Thus, with y 
equal to the ratio of specific heats, assumed constant throughout, the sonic speed a, 


at r = r, is given by 
a; = YP; / pr and Qrp,ayr, = mM. (1.6) 


The continuity equation then becomes 


p*wr* = 1. (1.7) 
Here w is always positive for sink flow. The equation of state is now 
p* = 6p*. (1.8) 


Eliminating p and p in Eq. (1.1) by using (1.3) and (1.4), and introducing the non- 
dimensional quantities, with » = log r* as the independent variable, we obtain 


att“ 
an Y Lan \W Ww 


2 . . * 
- ~o6*4 u*(1 + n( oe = v) + la + k) = + bo | du") (1.9) 


where 


2r 
a* =(R)' == S2, ip’ = (1+ 8h. (1.10) 
piQyr; m 
a* thus denotes the inverse of the Reynolds number and will be considered much smaller 
than unity throughout this paper; whereas k expresses the relation between the two 
viscosity coefficients. Stokes’ assumption on the value of yu’ states that 


p’ = 0 or k = -—1/3. (1.11) 
As this assumption does not agree with observations for many kinds of fluid [6], condition 
(1.11) will not be imposed on the final calculation of the flow field. 
By using again Eq. (1.3), the energy equation (1.2) can be integrated once to yield 
the following nondimensional form 
o dé 


we x 6 : * {<< # * dw* 9]. | ee + a 9 
= Pe teat yo ia t t+ : + 2ku = —D’ (1.12) 





where o denotes the Prandtl number 
o = Cyy/d. (1.13) 


The integration constant on the right hand side of Eq. (1.12) is chosen as shown above 
so that the limit solution for vanishing viscosity agrees at large r with that of a nonviscous 
iso-energetic flow. 

Equations (1.9) and (1.12) are the two equations for two unknowns w and @. The 
boundary conditions for them can be determined by requiring that they tend to their 
respective inviscid solutions as 7 © so that these two solutions can appropriately be 
compared later. 

1.1 The inviscid solution. The solution for sink flow of a compressible inviscid gas 
can be literally obtained by putting a* = 0 in the above equations without justifying 
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the validity of such a simplification. The solution of this reduced system of equations is 


known to be 





-" 1 ae 1 ¥ -1/(y—1) 
- =e’ =w (4 -— 1— wu) ’ (1.14) 
: 2 2 
and 
1 —l1, : ie ; 
g=2ti_1 sw, «=p = (p*)”" = Orr. (1.15) 


The value of r,; can be expressed in terms of stagnation state as 


m m / 1 (1/2) (y¥ +1)/(y-1) 
sce ~ 2dape ( 2 " : (1.16) 


27rd; p; 277A Po 2 


where a) = YPo/po and po , po are the isentropic stagnation pressure and density. Equa- 
tion (1.15) simply states the iso-energetic and isentropic relations. 
This inviscid solution w(r) given by Eq. (1.14) is plotted in Fig. 1. It gives no solution 


forr <r, ; but forr > 17, , wis a double-valued function of r. On one branch w tends to 


} +i ee 
{4-1 


























r=/7, 


Fic. 1. Graph of inviscid solution. 


zero so that thermodynamic variables tend to their stagnation values as r >; on the 
other branch w tends to the maximum speed attainable [(y + 1)/(y — 1)]'”, and the 
thermodynamic variables tend to zero as r >. They will be designated as subsonic 
and supersonic branch respectively. Both of these branches terminate at r = r, with 


sonic speed (at which the fluid speed equals the local speed of sound). The slope of the 
curve w(r), 

dw 9 hana aly + 1 y—-1 ,\“r-» 

er “(w = 1) — SS wo . 


= ——— w — 

drt vyt+1 2 2 
is much smaller than unity for r > r, on both branches and consequently viscous effects 
become comparatively unimportant there. But as r — r, , w — 1, dw/dr becomes nu- 
merically unbounded, and thus the viscous force near r = r, should play a role as sig- 


nificant as those of inertia and pressure forces. 
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Now this inviscid solution will be used as a guide to study the sink flow of a real 
fluid governed by Eqs. (1.9) and (1.12) for large values of Re, in the sense that it is 
assumed that the limit of the viscous solutions for vanishing viscosity approaches the 
inviscid solution as r +, for both subsonic and supersonic branches. By continuing 
these viscous solutions backward in r where viscous effects become more and more 
prominent, it is expected that the real fluid, affected by viscosity and heat conduction, 
will flow across this fictitious sonic circle, which is a limit line when viscosity is neglected. 

It may be remarked here that the equations for source flow of real fluid can be obtained 
from (1.9) and (1.12) by changing the sign of the terms with factor a* if w again represents 
the absolute value of the radial speed, normalized relative to a, [5]. Hence the inviscid 
solutions for source and sink flow are identical, but their respective viscous solutions 
will be shown later to have quite different features for all r. 

2. Properties of the solution curves. 

2.1. Approximate differential equation in the “phase space’. In order to study the 
qualitative properties of the solution curves, several assumptions will be introduced in 
this section to simplify the analysis while most of the important features of the original 
system will still be maintained. The Prandtl number, ¢, is assumed constant because u and 
\ have almost the same dependence on temperature. In this section, u is also taken to be 
constant so that u* = 1. When the complete solution is calculated later, these assump- 
tions introduced here become unnecessary. 

Equation (1.12) can be integrated when the Prandtl number 


o = (1 + 4a*k)/(1 +k) (2.1) 


(under Stokes’ assumption, k = —1/3, then ¢ = 3/4 — a*), and the final integral is 


6 — }? TS ico (1+ ta*i)u" = E exp (—on/a*) = x(t)” (2.2) 
2 2 : r; : . 
where FE is the integration constant. The value chosen above for ¢ is actually not far 
from experimental data (o = 0.72 for air at standard condition). As ¢ only appears in 
the coefficient of the derivative in Eq. (1.12), it follows from the theory of differential 
equations [Ref. 7, p. 142] that the solutions and all their derivatives will be continuous 
ino forw > 0, — © <o <o. Thus the assumption of choosing this particular value of 
o« would merely lead to simplification of analysis rather than any material change of the 
solutions. If we further require by physical argument that the deviation from the iso- 
energetic relation expressed by the term with the arbitrary constant F shall not over- 
whelm the left hand side terms for r < r, , we may assume that E = 0. This restriction, 
however, can again be relaxed when the complete solution is discussed later. It will then 
be shown that FZ is indeed of the order 0(a*). Thus the particular solution with E = 0 
would still provide a good approximation to the complete solution. 
Substituting Eq. (2.2) with EF = 0 into (1.9) and eliminating the explicit dependence 

on 7 by the substitution 

> dw dw 

V= - dn me sey’ a ’ (2.3) 


we obtain 


ie V 
dw 


aw’ } 


+ V[l — (1 — aa)w’] — w{l — [8 + (a — lajw*} = 0, (2.4) 
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where 
by > Sine = 3 ‘ 
= (] La*, 8 = - = ——— «o ye 
+1" yt1’? “14k ¥ ini 


The variable V is closely related to the fluid velocity gradient. Since the term aa and 
(a — 1)a@ in the brackets are merely corrections to constant coefficients of 0(1), the 
properties of Eq. (2.4) would not be altered if we had neglected these terms in order to 
simplify further algebra. Thus the approximate differential equation 


aw'y 2 + V(l — w’) — w(l — Bu’) = 0 (2.6) 


du 


in the phase space (w, V) is expected to exhibit all important features of the original 
system, Eqs. (1.9 — 1.12), for u* = 1. An equation similar to (2.6) was derived by 
Sakurai [4] and later was discussed in detail by Levey [5] for source flow in a real fluid. 

2.2. Properties of the solution curves in phase space. Equation (2.6) is nonlinear and 
cannot be integrated. However, several important features of the solutions can be readily 
seen by studying the properties of the vector field (w, V) defined by Eq. (2.6), such as 
the type of its singular points, the curves of zero slope and zero curvature together with 
some obvious isoclines. 

2.2a. The curve of zero slope; the inviscid solution. Let G, be the curve on which dV /dw 


given in Eq. (2.6) vanishes, G, is then given by 


V = w(l — Bu’)(1 — w’)", (2.7) 
which is also the inviscid solution in w — V plane. The function V,(w) given by Eq. 
2.7) has a simple pole at w = 1 (the fictitious sonic circle), and two zeros at w = 0 and 


w = 6” which correspond respectively to the subsonic and supersonic branch at 
r= o, Near the origin, V,(w) given by Eq. (2.7) has the following power series expansion 
= 2.8) 


VY, =w+(1— ~w+(1-— pw’ +(1-— Bw +-:--, 


; ; ; ; ; ; 1/2 : 
which starts from w = 0 with slope unity. Near the point w = 6 “’, the expansion of 


V,(w) is 
. 26x 1+ 38 2, , 1+ 68+ 8B ,21252 
v, = 8 [ —oq —-p% PT oq —p? Y ® 
1 + 108 + 56° , 12 43 - 
a 11 — 8! (8 x) -t | 2.9) 
where 
z=w—p'”, (2.9a) 


This branch of the curve G, starts from (8"'’’, 0) with the slope 28/(1 — 8). The curve 


G, divides the infinite strip 0 < w < 8” into regions of positive and negative slope as 
shown in Fig. 2. In this strip, dV/dw given in Eq. (2.6) becomes infinite only on (i) 


V =0,0 <w < 8’, and (ii) w = 0, V ¥ O. Besides, w = B''’, V ¥ Ois also an 


isocline on which 


dV/dw = (1 — B)/a. (2.10) 
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Fig. 2. Regions of positive and negative slopes; isoclines. 


2.2b." Properties of the singular points. The only singular points of Eq. (2.6) are (0, 0) 
and (6°'’*, 0). The origin is a singular point of higher order. But if one sketches the 
vector field (dw/dt, dV /dt) defined by 


dw/dt = aw’ V, dV/dt = w(1 — Bw’) — V(1 — w’) 


along a simple closed curve in the neighborhood of the origin with the origin in its interior, 
one finds that the Poincaré index (see [8], p. 45) of this singularity is equal to — 1. Thus 
the origin is a saddle point through which only two solution curves may pass. One of 
these is w = 0 which either yields a trivial solution (V = 0) or has no physical meaning 
(V # 0). The other solution curve starts from the origin with slope equal to 1 (which 
coincides there with the inviscid solution) and thus represents the only possible radial 
sink flow with stagnation at r = . By substitution of a power series into Eq. (2.6) 
(or by ordinary iteration), the asymptotic value, for small a, of this solution near the 
origin is found to be 


V=w+ [1 — 8) — alw*® + [(1 — B) — a(5 — 48) + 40°]w* 
(2.11) 
+ [(1 — 8) — a(14 — 168 + 36’) + 10a°(4 — 38) — 270*Jw’ + ---. 
The point (8°'*, 0) is a singular point of regular type. In the neighborhood of 
z=w-—s’” =0,V = 0, Eq. (2.6) becomes 
dV _ —26r + (1 — B)V + P(x, V) (2.12) 
dx aV + Q(z, V) 
in which P and Q vanish like x? + V’ as xz, V — 0. Thus the singularity [Ref. 8 pp. 


37-44] is (i) a nodal point if a < (1 — 8)’/(88), (for air, 8 = 1/6, a < 25/48), or (ii) a 
spiral point if a > (1 — 8)’/(88). As the problem will be confined to the case a < 1, 
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we shall only consider this singularity to be a nodal point (which changes to a saddle 
point for source flow [5]). All solution curves passing through this point will have at this 
point two distinct slopes which can be calculated from the secular equation of Eq. (2.12), 


namely, 


r —a 
| = 0, (2.13) 
28 A= Ck — 8) 
which has two unequal positive roots 
A, = (1 — B) — 2aB/(1 — B) + O(a’), A» = 2aB/(1 — B) + O(a’). 2.13a) 
From these two eigen-values two eigen-vectors associated with Eq. (2.12) at (x = 0, 
V = 0) can be obtained. By using this result it can be shown that the solution curves 
passing through (x = 0, V = 0) have near this point the following parametric repre- 
sentation 
28 1-6 2 
ia? = > A _ af 
V- x=Cee", J at — )e = Cee ; 2.14) 
1— 8 5 \ a l — £, 
Since \, > A. > 0, V, 2 —~O0ast— — o. It also follows from the above equations that 


there is an infinite number of solution curves, corresponding to arbitrary C, and C, 
(C, # 0) which have the asymptotic value 


8 
V= = 52 +Cl\z sat tog near (x = 0, V = 0); (2.15a) 
and, in addition, there is another solution curve (~ C, = 0) passing through this point, 
of the value , 
__ (1-8 28 | ; P 
V= | _ v near (x = 0, V = 0). (2.15b) 
a 1 — B/ 
The first group of solutions, given by Eq. (2.15a), have the same limiting value at r = © 
as the inviscid solution and hence represent the many possible sink flows starting with 
maximum velocity w = 6”? at r = &, while the solution given by Eq. (2.15b) is 
physically irrelevant. The asymptotic value, for small a, of this physically significant 
solution near the point x = w — B “3 = (is 
. 2Bx 2aB 8a°B° : 
VY & 4; 1-4 s+ ; + O(a’) 
l1—f8 1 — £B) (1 — B) 
Br’*x 8(3 + 138) ' 
_ (1 + 38) + 2a - s— + O(a’) 
2(1 — 8B) (I — fg) : 
(2.16) 
Bx” . 9 7 ee ap 2 , \ 
+ 3] (1 + 68 + B) + 2(5 + 426 + 3328) 3 + O(a’) | + ? 
2(1 — 8B) (1 — B) J 
—cisr 


where C is an arbitrary constant. The last term follows from Eq. (2.15a). 
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2.2c. The curve of zero curvature. The second derivative of the solution V(w) can be 
deduced from (2.6). Let G, be the curve on which d’V/dw’* vanishes, the equation for 


G, is, except where wV is zero, given by: 


2V* — w(l + Bu’)V*? + a '(1 — Bw’) [V1 — w’) — w(l — Bw*)] = 0. (2.17) 


The function V,(w) satisfying this equation has the following properties: 
(i) For small values of a, V, has only one real value for either 


_ 


- 3 ; fey 1/3 | Qa | 

w 4 /2 » —= =, * 
Sw<itgl B)a/2] or w>B 1+ iq BI 
V. has three real values for 


3 nii/3 ; -1/2 . a * 
+5 [0 — Ba/2l" <w <p E + 4(1 — 5} 





V. has two equal real values at 


(1 — B)a/2]'”; =p i+; + - 7 and w= ,”, 


w=1+ 


bo | Go 


(ii) The curve G, starts from the origin with slope equal to unity and has, in the 
neighborhood of the origin, the following expansion: 


V,=w+ [(1 — 8) — alw*® + [(1 — 8) — a(5 — 48) + 4a’ }w* 


(2.18) 
+ [(1 — 8) — a(14 — 168 + 36’) + 20°(17 — 128) — 21a*}w’ + 
(iii) The curve G, crosses w = 1 at 
V1) = (1 — 8)?°(2a)"'*[1 + O(a'”*)] (2.19a) 
with the slope 
dV./dw = (2/3)(1 — 8)'*(2a)"**[1 + 0(e’”’)]. (2.19b) 


(iv) When w = 6°'”, V; = 0, and V, = 6°”. Thus G, has a double point 


at (8°'*, 0), where the curve has two different slopes 26(1 — 8)~* [1 + O(a)] and 
~' [1 + O(a)]. Near this point these two branches of the curve have the following 





(1 — B)a 
expansions: 
—— - 208 , 8a°8 . 
V2 =i. A) [1+ — peta 28 + 0) | 
a 2 =| +38 ) + tn SEHD + oy 5] (2.20) 
21 — ai- 7h 
+ sq gr ll + 68 + 6°) + 0(@)] + +} 
and 
Oe det | _ 208 5] mo 9 
Vy” = —— wf 1 a — pt 0’) + (2.21) 
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(v) The slope of the curve G, becomes infinite at 
w~1+ ${(l — Ba/2]”", Vr a1 — 6/2)” 


and 


—1/2 — as Ro -1/2 
os ee aa 


P f 7 , 3 fc 
(vi) For w > 1, V; ~ Bu, V2 ~ Bw'/2. 
The curve G, divides the (w — V) plane into regions of positive and negative curva- 


tures as shown in Fig. 3. 





V | 
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Fic. 3. Regions of positive and negative curvature. 


2.2d. Sketch of solution curves V(w); definition of cylindrical shock. The above dis- 
cussions on isoclines, types of singularities, regions of positive and negative slope and 
curvature, enable the solution curves to be sketched. 

First let us consider the solution curve starting from the origin. Comparison of Eqs. 
(2.8), (2.11) and (2.18) shows that, for a and w both small, 


Vi(w) > V(w) > V.(w), (2.22) 


the difference (V; — V) being of 0(aw*) while (V — V2) ~ O0(a’w’). Thus the solution 
curve V(w) lies in between G, and G, where the slope and curvature of V(w) are both 
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positive. A careful study of the slope of V(w) and G, indicates that the solution curve 
lies above G, for w > 0. Hence V(w) is a monotonic increasing function of w with in- 
creasing slope, passing through w = 1 between the points V = (1 — 8)/a and V = 
(1 — p)?” (2a), and finally ending up at w = 6”'”” with the slope dV/dw = (1 — 8)/a 
(cf. Eq. 2.10 and Fig. 4). 

As previously shown in Eq. (2.16), there is an infinite number of solution curves 
starting from (w = 6”, V = 0) with the same slope 28/(1 — 8). However, 
for (w — 6B '”*) and @ both small enough, comparison of Eqs. (2.9), (2.16) and (2.20) 
again shows that all these solutions satisfy, for small negative (w — 6~'””), the following 
inequality 


Viw) < V3"(w) < V,(w), (2.23) 


as shown in Fig. 4. As w decreases from 6~'”’, V (for every finite C in Eq. 2.16) decreases 
with increasing slope until it intercepts G, with a positive slope. For further decrease 
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Fic. 4. Sketch of the solution curves in the phase space. 


in w, the curve V should lie above a straight line with this positive slope at the point of 
intersection because V has positive curvature in this region. Hence the solution curves 
will eventually meet G, with zero slope. From there on, for further decrease in w, V 
increases from negative values and later crosses V = 0 with infinite slope at some point 
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in between w = 0 and w = 6”, as can be shown by the method of bounding curves, 
and will be made explicit in our later calculation. Further extension of these solution 
curves shows that V increases with increasing w and finally approaches asymptotically 
to the subsonic branch solution which starts from the origin. There is a particular value 
of the integration constant C (in Eq. 2.16), say, C = C, < 0, for which the solution 
curve finally ends up at the origin with infinite slope. For C < Cy , the solution ceases to 
have physical meaning. On the other hand, the solution curves for C > C, have a very 
interesting feature that these viscous solutions all exhibit the transition process from 
the inviscid supersonic branch toward the viscous subsonic branch. Let us consider, in 
particular, the solution with C = 0. It first intercepts G, at (w, , V,), say, and then crosses 
V = Oat w = w,. Since w, > 1 and w, < 1 (as will be shown later), the flow between 
these two states may thus be defined as that of a “cylindrical shock”’.* Inspired by the 
result obtained in Eqs. (2.19a, b), we see that the equation governing such a cylindrical 
shock flow can be approximated by the following similarity transformation 


w=1+a'o, V =a’, (2.24) 
which will render all terms in Eq. (2.6) equally important in this flow region, that is, 
for vanishing a, 

d di/dw = 20 + (1 — 8). (2.25) 
Since this equation governs the flow of both branches near the sonic speed, it may be 


called the “equation for cylindrical transonic flow’’. 
2.3. Sketch of solution curves w(n) in physical space. From the definition of V(w) 


given by Eq. (2.3), we have 








n=— | [V(w)l* dw +h, (2.26) 
where the integral stands for an indefinite integral and 
} : | ( 2 (2.26 
ees og oa og | 2.26a) 


so that 7 tends to its inviscid solution for 7 large. 

Now for the subsonic branch, V(w) > 0, hence 7 is a monotonically decreasing 
function of w. Moreover, for some value of w, V(w) is less than its corresponding inviscid 
solution, V,(w) (see Eq. 2.22). Hence, from Eq. (2.26), 


Nvis(W) < Ninvis(W). (: .27) 
In other words, at every 7, (w),i, is slowed down from its inviscid value due to the 
viscous effect. 


*This terminology is adopted by both Sakurai [4] and Levey [5] to describe such type of flow. The 
term “shock’’ is borrowed from its conventional meaning to indicate the transition from one branch to 
the other, though the transition is rather different from that occurring in a plane normal shock. Perhaps 
this terminology relates more closely to the conventional meaning of a shock for the constant C slightly 
greater than C> (see Fig. 4), because then the jump in w and the slope dw/dy in transition become greater, 
and the position of transition is farther out from r = 7 (see Fig. 5). But since there is no adequate cri- 
terion to distinguish one from another value of C, we shall retain this name. Another terminology, the 
dissipation layer, is suggested by Prof. H. S. Tsien to avoid this ambiguity and, in addition, to stress the 
importance of viscous effects in this layer. 
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For the supersonic branch starting from w = 67’, V(w) < 0 for w, < w < B™”, 
hence in this interval 7 is a monotonically increasing function of w. Atw = w,, (~ 7 = 
n2 , say), (dw/dn) = — V(w.) = 0, and, from Eq. (2.6), (d’w/dn’) = (VdV/dw),.0 = 
(1 — Bw,)/(aw,) > 0, therefore w(n.) = w, is the only minimum of w on this branch. 
For » < m2, 7 decreases with increasing w. Furthermore, because of the relation given 
in Eq. (2.23), we have 


Nvis(W) > iavis(W) for w>1,V<0. (2.28) 


The above properties of the solution enable the solution curves to be sketched, as 
shown in Fig. 5. As the constant C decreases from zero, the minimum value decreases 
and the jump in w increases while the jump takes place farther upstream. 











y= beg 7 


Fig. 5. Sketch of the solution curves in the physical space. 


It follows from Eq. (2.24) that if we introduce 


n = at, w=1+a'd, (2.29) 
then Eq. (2.25) becomes 
a’ = d = 


which is the equation governing the flow near w = 1, 7 = 0. Integrating this equation 
once, we obtain 


dw/dt + w = (1 — BE + const. (2.31) 


This equation will be integrated and discussed in our final calculation. 
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3. Calculation of the solutions by using PL K-method*. In this section we shall 
calculate w(n), 6(n) governed by the original system of Eqs. (1.9-1.12). Throughout this 
section » will again be assumed constant so that u* = 1, but no restriction will be imposed 
on ¢ and k. Consequently Eqs. (1.9) and (1.12) become 


) ’ 32 
we o + E 6 - - ow =— i o( Fe - whe’, (3.1) 
dn  ¥ dy d 2y 7 








6+ 1 Sa (1 + baw? + ~ =| af +a-D ae | . a = 
2 4y > dn ay - 


where a is given by Eq. (2.5) and 
b= a/8B = (1+ 1/y)(1 + L/h". (3.3) 


Using the PLK-method, as described below, the generalization to the case un = u(T') 
presents no particular difficulty. The calculation for the case when yu is proportional to 
local temperature is carried out in Ref. 1, Sec. 7. The result shows that all important 
features of flow quantities for u4* = 1 remain up to the term of 0(a””’). 

It was seen before that, even for a simplified version of these equations, such as Eq. 
(2.6), the conventional perturbation method merely leads to asymptotic solutions for 
small a near w = 0 or w = 6” (see Eqs. 2.11, 2.16) because the coefficient of w* does 
not diminish as n —«. This asymptotic result fails to be a good approximation to the 
required solution as w deviates further from w = 0 or w = 8°” and becomes almost 
useless for calculations near w = 1. Now let us resort to the PKL-method, which is in 
essence to expand the solution in terms of power series in a with coefficients as unde- 
termined functions of a parameter é. 


wk) = — + aw" (—) + a’w(&) + +, 
né) = 98) + an?) tan (E) +, (3.4) 
A) = OE) + ad) +o OPE) + +e. 


The need of a parameter é to represent the solution and that w starts with the term & 
are clearly suggested by our previous discussions. Substituting these expansions into 
Eqs. (3.1) and (3.2), noting that 


dw/dn = w’/n’, d6/dn = 0’/1’, d’w/dn? = w''(n’)* — w'n’'(n’)’, (3.5) 


where prime stands for d/dt, and then equating equal powers of a, we obtain the zeroth 
order equations as follows: 


E +2 ‘ (0° /t) — £6 9” ey ”’? = 0, (5.6) 


g ee (xt I ae el f) |, We oe 0. (3.6b) 


If we choose 7°’’ different from zero, then we have the zeroth order solution 


g(g) eS a ion I £*, (3.7a) 


| | | 
n'(G) = —logé - a5 log | MEF - WS. acid 


*See the footnote in the Introduction. 
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The integration constant in Eq. (3.7b) has been so chosen that when a = 0, = w and 
n°’ (w) agrees with inviscid solution. 
The coefficient of a in the expanded equations gives the following first order equations: 


1t% 


ge? + (y — lew” = — 1—lip + —— (1 — Be’)(1 — 7s (3.8a) 


(nn + ow + (78? — 0) w” + ede — 505 
(3.8b) 


a (eo? 4 9) g"h ~ 7 + L ry (0)’ a £(n” ) l, 


where b is given by Hq. (3.3) and 


= ik | eo oe 
e= ; E 2e(1 + 5|- (3.8¢e) 


Substituting Eqs. (3.7) and (3.8a) into (3.8b), and rearranging the terms, we finally 
obtain [see Ref. 1, p. 25]: 


sin w ok ae 
(3.9) 
26 


pay 
aw ~ 0° 


+B 





_ 2 vé | 
is latet et +4 Gay 
where 


= S ee me = ” _ (1+ 628 
A= og? ( a)(1 B), B aI 3 (a — «) ‘rs ¥ | (3.9a) 


and the constants a, b, € are given by Eqs. (3.3) and (3.8c). Now in order to ascribe to 
w" a possibly lowest order singularity at = 1 to improve the convergence of the series, 
we decompose Eq. (3.9) as follows: 














dé [nw] = A ane Dp? + B af Dp (3.10a) 
‘ = -0 — oe ~— 0° + 0 - 0"). (3.10b) 

The solution of Eqs. (3.10a, b) is then 
w°(g) = — Ag(1 — &)"" — Bel — Be)(1 — &°)* log | us —— (3.11a) 
7° = -(1- ala — ef)" + : log | 1 — & 7 (3.11b) 
The integration constant of Eq. (3.9) is first absorbed in n‘” and is then omitted because 
of its negligible contribution near — = 1. In the above first order solution, w‘” and 


n have singularities at £ = 1 of the same order. After w” is so determined, 6‘ (£) is 
then given by Eq. (3.8a). 

Proceeding in a similar manner to obtain the second order equations by equating 
terms with a’, we find that the resulting equations possess solution of quite lengthy ex- 
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pression, in which »“(é) starts with the term 2(1 — 8) (1 — «) (1 — ¢’)~* followed by 
terms of 0((1 — £’)~*), while w‘ (é) still can be made to be of 0((1 — *)~'). Therefore 


the final solution can be expressed parametrically as follows: 


E ‘1 — Bt’) 1 —1., 
wt) = & — of. { 1 = e2 T ph log i —_— y— a | 
. 4 fe i 3. 12a) 
a 
+ (3) 
l et w#w-~ is 
n(é) = -| tore + — j log — — =¢ ] 
9 ane — 
—a(l — a| te + Flog l —F |+e te? 
, é (3.12b) 
( ‘ 3 ¥ ae | (1-8) */(2a8 
eT : “ss —_—_—_—_—— = 2 
+ 1B aes) + 2 = * 
where A, B are constants given in Eqs. (3.9a) and C is an arbitrary constant for the 
solution starting from ¢ = 6”, but C = 0 for the solution starting from & = 0. The 
term with C is included in the above equation in order to exhibit the behavior of solution 
curves near the nodal point (w = 6-'””, 7 = ), as previously discussed in Sec. 2.2b. 


This term is of higher order than any finite power of a. In fact, it has an essential singu- 
larity at a = 0 and hence cannot be obtained by the expansion procedure. The value 


of 6 is given by 


ag = (Vt1_ r= 1p) 


2 a 





P 1 . 1 —1., 

_— pet att a te tS ~ & (3.12c) 
_ by», « Hl — BE a = 6) | (t=) 
5 § + — ? +0 (1 — ¢’)’ . 


Several interesting features of the above solution may be mentioned here. (i) « = 0 
foro = [2(1 + k)]”’ (= 3/4 if k = — 1/3, see Eq. 3.8c). With the value of o and k 
lying in the experimental range, « is still a small number. Consequently, the variation 
in ¢ and k only contribute a small correction to the coefficients of 0(1) in the solution. 
This fact confirms our previous statement in Sec. 2.1. (ii) By substituting Eq. (3.12) 
into Eq. (2.2), it can be found that the arbitrary constant E is of 0(ae). (iii) The most 
important property of the above solution is that it does not provide an approximate 
solution with uniform accuracy over the interval of — such that 0 < w < 8°’. Asé 
approaches unity, the higher order terms, especially in n(¢), become more important 
relative to the zeroth order term. More — isely, the solution is good only for : <3< 


£ 


1 — Ka” and 1 + Ka’? < ¢ < Bp”, K being a positive constant of 0(1). At é = 1 
+ Ka'”*, all terms in the expression for (¢) become of the same order, 0(a*”*); red the 
convergence can be made sufficiently rapid by an appropriate choice of the value K. 

In the subsequent calculation of the solution through the transonic region, we shall 
only consider a particular solution with C = 0 and e = 0 in Eqs. (3.12). Furthermore, it 
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has been found convenient to take K = 2(y + 1)~’”*. With this value of K, we obtain, 
from Eqs. (3.12), the following result: 
(i) on supersonic branch, at w = 1 + 2a'” (y + 1)", 
n = 2.28(y7 + 1)'%a"”" and dw/dn = 0.615a7~'"(y + 1)~*”; (3.13a) 
(ii) on subsonic branch, at w = 1 — 2a” (y + 1)°™, 
n = 1.766(y7 + 1)'“a’*” and dw/dyn = —0477a'"(y + 17°”. (3.13b) 


These values will serve for the boundary conditions imposed on the transonic solution 
to be obtained below. 

That the PLK-method is powerful in solving this problem can still be stressed further 
by the following argument. As the first order term in the expression for w() (see 3.12a) 
is quite unimportant in the aforementioned regions of ¢, one perhaps would try, instead 
of Eq. (3.4), a simpler expansion 


n(w) = 9°°(w) + an’?(w) + --- (3.14) 


and a similar expansion for @ in terms of w. It can be shown that the above expansion 
will yield a solution in which n°’ is identical to inviscid solution, but 7” has, in addition 
to a simple pole at w = 1, a pole and a logarithmic singularity at w = 6~'””. Consequently 
the assumed expansion (3.14) becomes invalid for r large on the supersonic branch, and 
thus leads to an erroneous result. 

4. The solution in the transonic region. To obtain an approximate solution in the 
transonic region, as discussed in Secs. 2.2d, 2.3 (see Eqs. 2.24, 2.25 and 2.29) and also 
as guided by the boundary conditions (3.13), first we distort the independent variable by 


the transformation 





n = at (4.1) 
and then expand w, @ into the form: 
w(t) = 1 + a? w() + a’? w () + aw (&) + s+, (4.2) 
af) = 1+ a? OME) + a0 (6) + aed) + +e (4.3) 
With « = 0 (~ 2e(1 + k) = 1), Eqs. (3.1) and (3.2) reduce to 
. 2 
i a (hn aa)w?| + oll — 64 — Deel oa’ SS, (4.4) 
” dn 
tot + 9 3 LL + baw? + O(a’). (4.5) 


Substituting (4.1) and (4.2) into (4.4), we obtain the first order equation: 


2 (1) (1) 
a 4+ 2y “ =(i-p (4.6) 
Ss 


and the second order equation: 
dw” d ' d 
av = toy. & (1)y3 (1) 5 
de +2 dé (w'’w"’) dé (w*’’) (1+ Bw’. (4.7) 


The correction due to the terms with constants a and b enters only in w® 


and higher 


order terms. 
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Now Eq. (4.6) can be integrated once to yield: 


- y =x+n, (4.8) 
dx 
where 
1/3 —1/3 
, , l : 
y = (r+ ‘) ., r= (r+) f, (4.9) 


and the integration constant x, can be determined by using the boundary condition 
(3.13) to give 
x, = 0.033 for the supersonic branch, (4.10) 


x, = — 0.0056 for the subsonic branch. 


The nonlinear equation (4.8) (see also Eq. (2.31)) is an approximation to the Navier- 
Stokes equation, a case in which all the inertia, pressure and viscous forces are equally 
important. Now Eq. (4.8) is of Riccati type, which, by using the transformation 
1 dv . 
yz) = -—, (4.11) 
v dx 


can be reduced to a second order linear equation: 
d’v 
5 — («x +2, = 0. (4.12) 
dx 
The solution of this equation for (x + x,) > Ois 
v(z) = Me2'“[I_1/3(2) + NI,,3(2)], 2 = 3x + 72,)°”’, (4.13) 
where J(z) is the modified Bessel function of the first kind, and M, N are the integration 


constants. By using Eq. (4.11), the solution of Eq. (4.8) for (x + x,) > 0 is then 


oe Ni ale) : are 
x) Toj3(@) + NI-2/s(2) z==(¢+2,)*” (4.14) 


y(z) = (# I wld + NI, 3(2) ' me ’ 


Col bo 


The constant N can be determined by using again the condition (3.13). 
The continuation of Eq. (4.14) into the region (x + x,) < 0 is provided by 
2°. n@ = —E" IAs), 2 Tan) = "I an(d), — 
(4.15) 
27, (2) _ 7 Jo(0), gly ,(2) _ ae 


where 
¢ = 2°"? = 2[-(r + 2,)]*”. 


Consequently Eq. (4.14) becomes, for (x + 2,) < 0, 
3t\'" Jo,3(¢) + NJ-2;3(8) 2 a 

yo) = (4) Re t= -[(-(¢ + 2]. (4.16) 

ms 2 J ~1/3(9) iit NJ1/3(0) . 3 | ( * JI 
To discuss the above solution, we first note that the inviscid solution in this transonic 
region, to the first order approximation within the O(a’”’), is 
2 = 2, (4.17) 
which has two branches for x > 0 and gives no solution for x < 0. Now before we de- 
termine the value of N for the corresponding viscous solution, we may also note that the 
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general solution, given by (4.14) and (4.16), is a semi-transcendental function of x, and 
the second integration constant N. It can be shown, from the property of J,(z) at large z, 
that in Eq. (4.14) 


youn” as r—@ when N > -1, 
y——2'” as «© when N = -l, 


and y has a simple pole at a certain finite z for N < — 1 (which is, of course, of no physical 
significance). This result shows that the viscous solutions tend to their respective inviscid 
values for x large in a manner which implies again that (w = 6°’, » = @) is a nodal 
point (admitting more than one value of N) while (w = 0, 7 = @) is a saddle point 
(admitting only one value of N). However, for zx + x, < 0, Eq. (4.16) shows obviously 
that y(¢) has an infinite number of isolated simple poles at ¢ = ¢, where the denominator 
vanishes. Since the properties of the solution curves in the (w, V) phase space exhibit 
no such singularities, the solution (4.16), therefore, represents a good approximation to 
the real flow only for ¢ lying in the interval 0 < ¢ < ¢, — 6 < ¢, , where ¢, is the first 
pole and 4 is a positive number, appropriately chosen such that y(¢,; — 5) is not yet too 
large to void our approximation (4.2). 

Having obtained the first order solution w“? (¢) = 2'”° (y + 1)~’” y, the second order 
equation (4.7) can be then integrated to yield 


w(t) = ee il wie? dé, (4.18) 
where 
g(t) = [ w'() dé and v(t) = [w’®]* — (1 + Bee) + const. 


It is obvious that w() is bounded wherever w“?(£) is bounded. Consequently, the 
approximation is good even if we only take the first two terms in (4.2) and (4.3). 

In order to obtain some numerical results, we first determine the value N in Eqs. 
(4.14), (4.16) by using conditions (3.13) to obtain 


N = — 0.585 for the supersonic branch, (4.19) 


N=-1 for the subsonic branch. 


With these values of x, and N (see 4.10 and 4.19), the solutions are plotted in Fig. 6 
(by using tables, Ref. [8]) from which several interesting results can be deduced as 


follows: 
(i) For the supersonic branch, the transonic solution starts from point A (see Fig. 


6) with the coordinates 
w, = 14+ 1.6[(y + 1/2) '"%a"™, m = 2.9[(y + 1)/2]'%a’”. (4.20) 
After the solution curve passes through a point of inflection G and then crosses the line 
y = 0 (w = 1) at point B (see also Fig. 4) where 
2/3 _ —1/3 


ne = 1.02[(y + 1)/2]"7a"(~V_ = —1.05[(y + 1)/2]°7a"™, 
(dV /dw) » = —0.95[(y + 1)/2)'"a"*”*), 
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Fic. 6. The flow velocity in the transonic region. 


it reaches a minimum when it intercepts the curve y* = x at C where 
w. = 1 — 0.45[(y + 1)/2)7'a'”, no = 0.2[(y + 1)/2]'a’”. (4.21) 


It then increases from w = w, to w = 1 at point D where 


Ip = ~o.sa(7+ 3) * 


(~v, = = (). o(%t le 1)" on”, dV ‘fn 1.1(r+ ) a), 
dw 2 


That is, there is an expansion wave following the cylindrical shock. 
(ii) The thickness of the cylindrical shock, as defined in Sec. 2.2d, is 


1\" as , 
An = ™ — 2 = 2.(>+2) es (in r-space, Ar = r,An), (4.22) 


across which the velocity has a total variation 
Aw = w, — w, = 2.05[(y + 1)/2)'"a"; (4.23) 
moreover, 
ww, = 1 + 1.15[(y + 1)/2)70"™. (4.24) 
Combining Eqs. (4.22) and (4.23), we obtain 
(Aw)(An) = 5.5a. (4.25) 


Comparing these results with those of a plane normal shock (e.g. Ref. [10]), we note 
first that the plane shock strength (~ Aw) is quite arbitrary while for a cylindrical shock, 
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Aw ~ 0(a'”*). The expression for shock thickness (Eq. 4.22) shows that An ~ O(a”), 
although, after combining with Aw, the expression (4.25) agrees with that of a plane 
shock [10] within the order of magnitude. The result (4.25) differs, however, from Levey’s 
result for the diffuse shock in a source flow (Ref. [5], Eq. 4.9), in which he explains the 
discrepancy as due to some degree of choice of the definition of the shock thickness. 
Our result also indicates that the maximum velocity gradient inside a cylindrical shock 
is of order a '”*, (in contrast to Levey’s result: 0(a~*)), while for a plane shock, the 
maximum gradient is of order (Aw)* a~' [10], which reduces to 0(a~'”*) if Aw ~ 0(a™*). 
The expression (4.24) differs from the Prandtl relation of a plane shock by a term of 
0(a'’*) which here agrees with Levey’s result (Ref. [5], Eq. 4.13). The differences between 
the present results of cylindrical flow and those of one dimensional plane shock can 
perhaps be realized by visualizing that the viscous forces exerting on the surfaces rdé 
and dr of a fluid element are indeed of quite different nature from those exerting in plane 
shock flow, since in the former case, the normal stress acting on the surface dr will have 
a component in the radial direction. 

(iii) On the subsonic branch, w is a monotonically decreasing function of 7. w = 1 
at point FE (see Figs. 6 and 4) where ne = — [(y + 1)/2]'” a*” and the velocity gradient 
(dw/dn)e = — Ve = — [((v + 1)/2)'*" a’ (see also Eq. 2.19a), which shows that the 
solution curve of the subsonic branch passes w = 1 (see Fig. 4) slightly above the curve 
G.. 

(iv) The thermodynamic variables in this flow region can be deduced from Eqs. 
4) 


(4.5), (1.7) and (1.8). That is, in the expansion (4.3) and 
p*(é) = 1 + ay (£) + a’ (€) + eee 
(4.26) 
pM) = 1 +a p'P@ + a9 + ---, 
we have 
0@ = —@ — n(*44) ve, (4.27) 
p? ia yp? = Y g?, (4.28) 


et 


where y(£) is given in Eqs. (4.14) and (4.16). The value of p‘” and 6” vs. 7 is plotted in 
Fig. 7. The supersonic branch starts with compression and is then followed by an ex- 
pansion wave, while the subsonic branch expands continuously. Equation (4.28) simply 
states that p*, p* and @ satisfy the isentropic relation up to 0(a'”*). This implies that the 
entropy variation, if any, across this region must be of order at least a’”*. 

It should be pointed out here that Eq. (4.17) fails to be a good approximation for 
n > O(a’). Consequently the transonic solution, which approaches the asymptote 
given by Eq. (4.17), cannot be extended beyond the transonic region. It then becomes 
obvious that the patching of these solutions with the outside solutions must be made 
at » = O(a’). This behavior of the present solution is further indicated by the above 
result that the values of the variables given by the outside expansion at 7 = O(a”) 
fit right into the transonic similarity rule. 

5. The solution in the inner supersonic region. The above transonic solution 
shows that the supersonic branch flow approaches asymptotically to the subsonic branch 








T. YAO-TSU WU [Vol. XIIT, No. 4 









































414 
2.0 a T wl T 
F-(e) 9", o-) &” Z cal 
- ~_ 
-T ail 
a ai 
= 
1.0 ponies 
| 7 & sussonic 
f sussomic —| _ 
CO 
ZS s x= (ieay 
r) 
~~] |- § supersonic 
| 


_ 


F | | \ hak 
P supersonic _ | 


(a: lees gama Dames "i 
F ai aa | 


$+ cae i] —_ ie 
| lig, 


== 


| 
-2.0 - 1.0 




















+4 




















1.0 2.0 3.0 


° 


Fic. 7. The thermodynamic variables in the transonic region. 


and for 7 < nz both branches have supersonic local speed. We shall proceed to find the 
solution for 7 < nz . Let us consider first the continuation of the subsonic branch. For 
the sake of convenience, we shall take the point EF in Fig. 6 as the boundary such that 


at E, wz = 1 and 
“ y¥+1)' ” wr , . y¥+1 . 1/3 dV ° ++ 5 oie 2/3 ( 
"z= re] a’, Ve = "5 a , i = — & ‘ (5.1) 
a “ aAW/ sz “ 


Now rewrite Eq. (2.6) in the following form 


dv _1 (1 = ) , l= bl, (5.2) 
dw a Ww aw J 


From the properties of the solution of this equation (see Sec. 2.2) we recall that for 
n <nz,V > Vz, that is, in this region V is at least of 0(a~'’*). It then seems convenient 
to adopt an iteration method to approximate the solution for 7 < n¢ . As the first itera- 
tion, we neglect the second term on the right hand side of (5.2) and integrate, using 
(5.1), to obtain 


VY = (w — 1)" + (1+ ye. (5.3) 


aw 2 
Substituting this value of V into the right hand side of Eq. (5.2) and integrating twice, 
using again Eq. (5.1), we obtain, as the second iteration, the following asymptotic result, 
for a small, valid for w > 1, 


1/3 ee, 1/3 7 
= 02 — (xt t) a’’* arctan | Be 2 (ut 1) a | + O(a log a). (5.4) 


~ 
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This solution is in good agreement with the value given by Eq. (4.16) for nz < 1 < ne. 
It can be seen at once that as w — 8” (the maximum velocity at which 6 = 0), 
approaches its smallest value 7,, , say, where 


Mm — Ne = —1.6[(y + 1)/2]'a"”. (5.5) 


Thus we see that the flow supposedly terminates itself at a distance of 0(a°”*) to the 
inner side of » = 0, beyond which there is no solution to our present system of equations. 
If one were to investigate further the possibility that one could still obtain a solution 
of physical reality for 7 < 7, , one would face some rather dubious situations. For 
instance, near 7 = 7,, , the density, temperature and pressure all become so low that 
the validity of the equation of state for a perfect gas (1.4) is questionable. Besides, the 
fact that the viscous stresses reach the magnitude of the fluid pressure near yn = 7,, sets 
a likely limit on the applicability of the Navier-Stokes equation (1.1) and also raises a 
question as to whether Burnett’s higher viscous terms (Ref. [11], p. 271) should be em- 
ployed to overcome the present difficulty. Of course, it would seem plausible to continue 
our solution further inward by assigning appropriate values to the arbitrary constant 
C in (2.16). Nevertheless, it is still impossible to bring the flow to » = — © (r = 0) on 
account of the singularity that pu ~ r-' near r = 0 (see Eq. 1.3). To clarify these rather 
vague points is beyond the scope of this paper, although such clarification is certainly 
desirable. 
6. The entropy variation. We define S to be the specific entropy, 


T dS = CpdT — p'' dp (6.1) 
then the energy equation (1.2) can be written as 


S ’ 
puT 7 = div (A grad T) + ®, (6.2) 


where ® is the viscous dissipation function, which in this case is, 


eae ly u\” 4,,_ ,udu- 
nies 3 sili 2 (4) + (“) | * 3 (u u) r dr (6.3) 


The above definition of S of a fluid element is clearly for an open system since the heat 
exchange by conduction, and hence a net flow of entropy, occurs with the neighboring 
elements. Thus Eq. (6.2) merely expresses the energy balance, in terms of S, of a fluid 
element—a system not isolated, in the thermodynamic sense, from its surroundings. 
The analysis of formulating the second law of thermodynamics for the fluid flow case by 
making the system closed has been investigated in some detail by Tolman and Fine 
[12] and discussed later by Curtiss and Hirschfelder [13] from the point of view of 
statistical mechanics. Their idea is, in essence, to state that the change AS in the entropy 
of a system should consist not only of the net increase in entropy produced by irreversible 
processes taking place inside the system, but also of the entropy carried into the system, 
due to conduction of heat energy, equal to div [A7"~' grad 7] per unit volume. Following 
Tolman’s notation, we may thus write 


DS DS > , 
Aor) 7 Ai ).. oF ( on r). (6.4) 
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In the present case, Eq. (6.4), after combining with (6.2), becomes 





. ds r P no ’ 
uT( = — (grad T)° + ©. (6.5) 
p dr irr. r ss : 

As we are only interested in the qualitative features of our later results, we simplify 

these equations by using the assumptions: 


, = |, ‘= 0, C, , C, = constant, o = 3/4. (6.6) 


= 2 ‘yy 


Then the nondimensional form of Eqs. (6.2) and (6.5) are respectively 


_6 ds_4y | _@6 Saal dws a 
 - | dn = 3 aes y— 1 dn’ + dy u dy + Ww (6.7) 
i), PO +E) -eee} a 
: aed l bes irr iu 3 ‘ 7 = | dn i dn = dn so w (0.5) 
where 
s = afc. (6.9) 


Though the sign of the terms on the right hand side of (6.7) is in general indefinite, the 
value of the right hand side terms of (6.8) is, however, positive definite. Therefore (s) ;,, 
increases monotonically along the fluid flow, as predicted by the second law for a closed 
system. Subtracting (6.8) from (6.7), we obtain an equation which can be integrated 


to yield 
ie wee ee i ig (6.10) 
3 dn 


where the constant of integration is so chosen that both s and s;,, tend to s as 7 ~~, 
8 being arbitrary. 


In order to see that s of the shock type flow reaches a maximum near w = I, we 
substitute Eq. (2.2) with E = 0 into (6.7) and obtain 
6 ds 4 dw lw 
——_ ~ = q+ (f2 = — v). (6.11) 
y—ldyn 3 dyn dn 


This equation shows that for n outside the transonic region, the variation in s is at most 
of 0(a). Within the transonic region, d’w/dn’, being of 0(a~'), overwhelms the rest of 
the terms in the bracket and hence (6.11) reduces to 


6 ds _ 47 4, ¢¥ ans a 
5 ay = tO Gg (1 + 007%). (6.12) 





It then follows that s assumes its maximum value at the point where the curvature of 
the w = w(n) curve vanishes (d’w/dn’ = 0 at point G in Fig. 6 and at this point d’s/dn° 
is less than zero). However, from (6.10), the quantity s + (4 ya*/3)d log 6/dn does not 
have an extremum in the entire flow region. The above result is very much the same as 
that of a plane shock, as it can be shown that the velocity has a point of inflection at 


sonic speed where the entropy is also a maximum. 
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Integrating Eq. (6.11) with the aid of Eqs. (3.1) and (3.2) under condition (6.6), we 
obtain 
s — & = log [0(wr*)‘'’~”] = log 6 + (vy — I(log w + 7»), (6.13) 
where 
So = log (popo”), 
so that s > s, as r — ©. This equation is actually the definition of s usually given for 
a perfect gas. Substitution of the solution (3.12) into (6.13) shows that 
As ~ O(a) for n > O(1). (6.14) 
Within and around the transonic region, we substitute the solution (4.1)-(4.3) into (6.13) 
and simplify the expansion, then we obtain 
8 — 8 =a (y — I)[(y + 1)/2]'"(a — y’) + 0), (6.15) 
where x and y are defined by Eq. (4.9). The value of y = y(x) is given by Eqs. (4.14), 
(4.16) and also plotted in Fig. 6. Equation (6.15) is consistent with the fact that s = s, = 
constant along the inviscid solution y* = x. The variation in s along supersonic and 
subsonic branches of our solution follows directly from the data shown in Fig. 6. The 
result is plotted in Fig. 8. As » decreases along the supersonic branch, the entropy s 
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first increases until it reaches the maximum s, + 1.2 a’ (y — 1) [(y + 1)/2]’”* at 
point G, then decreases and later assumes once again the value s, (the value of s at 
” = ©) at point C where w'” is minimum. After that s decreases rapidly with further 


decrease in » and eventually tend to — ~ as the flow solution terminates. On the subsonic 
branch, s decreases monotonically with decreasing 7. However, by substituting Eqs. 
(6.15) and (4.3) into (6.10), it can easily be shown that (s;,,.) increases monotonically 
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with decreasing » and the variation in (s;,,.) is of order O(a). Consequently, the result 
that s— — © as7— nmin can be explained by visualizing from Eq. (6.10) that d(log 6) /dn 
decreases beyond all bounds as 7 — nui, . Physically, this probably implies that the flow 
is rather far from its equilibrium condition due to the large velocity gradient, inducing 
a rapid decrease in temperature, which even the important heat conduction in this 


region cannot compensate. 
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ON THE STABILITY OF THE SPHERICAL SHAPE OF A VAPOR CAVITY 
IN A LIQUID* 
BY 
M. S. PLESSET anno T. P. MITCHELL 


California Institute of Technology 


1. Introduction. It has been shown by G. I. Taylor [1] that a plane interface be- 
tween two fluids of different densities in accelerated motion is stable or unstable ac- 
cording as the acceleration is directed from the heavier to the lighter fluid, or conversely. 
This stability analysis is limited to small amplitude perturbations of the plane interface; 
and it is found that a small distortion of the interface begins to grow exponentially 
with time in the unstable situation and to decrease exponentially in the stable situation. 
While experimental observations agree well with the theory in the small amplitude 
limit for which the theory is valid, it is known that there are significant deviations in 
the rate of growth of distortions in the unstable case when their amplitudes become 
appreciable [2]. 

It is of interest to consider the analogous stability problem for the case of a spherical 
interface between two immiscible fluids of different densities in accelerated motion. 
For perturbations in the spherical interface of small amplitude, it may be shown [3] 
that the stability criterion deduced by Taylor for the plane case is subject to important 
modifications. 

The stability problem in the spherical case may be formulated as follows. A fluid of 
density p; is contained within a sphere of radius R; a fluid of density p, occupies the 
region exterior to this sphere. The fluids are supposed to be immiscible, incompressible 
and nonviscous. The equation of motion for the interface radius as a function of time, 
R(t), is readily determined under the assumption that the initial and boundary con- 
ditions are spherically symmetric. If the interface is distorted from the surface of a 
sphere of radius 2 to a surface with radius vector of magnitude r, , then one may write 


r,=R+ p Se (1) 


where Y,, is a spherical harmonic of degree n and the a,’s are functions of the time to be 
determined. The stability of the spherical interface may be established by considering 
whether interface distortions of small amplitude grow or diminish. More precisely, it 
is assumed that 

la.) «RO, 


and that terms of order higher than the first in a, and da,/dt are negligible. In such a 
linearized perturbation theory, the a,’s are independent of each other, and it may be 
shown [3] that they satisfy the following differential equation 


d’a 3 dR da 

= —_ J = 9 
dt R dt dt Aa, = 0 (2) 
with 


, — Inn = Ipo — (n + In + 2)piJd*R/dte — (n — In(n + 1)(n + 2)o/RE (3) 
—_ [nps + (n + 1)p,JR ’ ‘ 


where oa is the surface tension constant. 


*Received October 11, 1954. This study was supported by the U. 8. Office of Naval Research. 
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The spherical stability problem has application to the behavior of growing or collaps- 
ing gas bubbles in a liquid. Penney and Price [4] have carried out a numerical solution 
of the stability equation (2) for n = 2 for the case of a pulsating gas bubble in water 
with an internal pressure, p,; , in the bubble given by 


p,R°” = const. 


and with a constant pressure, po , in the liquid at a distance from the bubble. In their 
computations surface tension is neglected. The numerical solution showed that the 
distortion amplitude a, is much larger when the bubble is near its minimum radius 
than elsewhere. The problem to be considered here is a cavity for which the internal 
pressure, p; , is constant, in a liquid at constant external pressure py . These are approxi- 
mately the conditions in a vapor cavitation bubble in a liquid. An analytic solution for 
the stability equation may be found under these conditions. 

2. Solution of the stability problem. For a vapor cavity in a liquid, the vapor density 
p, may be neglected in comparison with the liquid density p, . The quantity A of Eq. 


(3) then becomes 


(n—1)dR o 
A = 5 - (n - I 1)( 2), (4) 
R dE n n+ n + 2) ok 
where p = pz» is the liquid density. The equation of motion of the undisturbed interface 
[3] is 
I dR 4 $ (az) =" i 2a Re (5) 
dt 2 \dt p 


Equation (5) may also be written as 


d| p3:(dR\| _ 4p24R | pi — po — 20/R 
dt E (4 y'] = 28 at | p 


which integrates, when p, — po is a constant, to give 
,(dR\? (aie) for _ : : 
on (: = 2 wietantied' a ( 93 3 2 32 a) 
ane] R ¥7 ) Ro dt 3 (Pi — p)(R° — Ro) — 4ro(R° — Ro), (6 


where R, is the cavity radius and dR,/dt is its radial velocity at t = ¢, . This integral 
of Eq. (5) is to be recognized as the energy integral of the system. 

The general features of the asymptotic behavior of the distortion amplitude a, may 
be made evident in a straightforward way. For the case of an expanding bubble, Eq. 


(sey 2P eo oz 
dt 3p’ , 4) 


with P = p; — po > 0, and it then follows from Eqs. (5) and (2) that 


(6) gives 


a, — const., R—- oo, (8) 


For the case of a collapsing bubble, it is convenient to transform Eq. (2) by the sub- 


a, = (Bo) (9) 
»=\p im 9) 


stitution 
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into 


d’b,, 
dt? 


' 3d|1dR 9/1 dR/’ 
OO” oa E dt | +4 E aR] +A 





— Gib, = 0 (10) 


with 


31, (R) O41 ER oy _ — 
= 7B (d + R de (n 1)(n + 1)(n + 2) ok (11) 
From Eq. (6), one has 
dkR|?  R: | (ate) 2p 22 | | | 
— | = — +. = D) 
E | R?\\at) + 3p 7 prj * LR (12) 
where p = po — p; > O for this case. The radial acceleration, d*R/dt’, is determined 
by Eq. (5); and the function G(¢) is found to be 
| 3n R3 (sy 2p 2e | 
QA) ~ —> pills = t 
O~ ~ oR | ai *s* mp = 7* aad 
except for smaller terms. It is evident that 
, me” 
G(t) R®? (14) 


where c is a real constant. One may now write a W.K.B. approximation to the solution 
of Eq. (10) for small RF in the form 


b, ~ [G(t)]"'* exp {+f (a(t’)]'” aw} ~ R** exp {sien [ Bag wh. (15) 
The distortion amplitude a, is then given by 
a, ~ R-™* exp {tion | Ro? a}, R—0; (16) 


so that a, increases like R~'* in amplitude and oscillates with increasing frequency as 
R — 0. This behavior has been found by Birkhoff [5] by a different procedure. It is of 
interest that the instability found by Birkhoff near R = 0 is qualitatively unaffected 
by surface tension. 

The question remains over what range of RF is the linearized perturbation theory 
for the distortion amplitude’, a, valid or consistent. The following problem will therefore 
be solved. A spherical cavity with radius R, at t = 0 expands, or collapses, from rest, 
dR,/dt = 0, under a constant pressure difference; at ¢ = 0, the cavity is supposed to 
have a distortion of small amplitude a, , and the subsequent behavior of a for any R 
is to be determined. Complete solutions for this problem are readily found when surface 
tension is neglected and these sqlutions are given first. The effects of surface tension will 
then be illustrated by some special solutions. 

(i) Expanding cavity; no surface tension 


*The subscript n for the distortion amplitude will be omitted in the following. 
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With no surface tension, the stability equation to be solved simplifies to 


da 3 dR da (n—1)dR 
— —— — - -a = 0; 7 
dé * R at dt R ae? % (17) 
one finds from Eq. (6) that 
dR | 2P Ro "= 
= - l1— gs 8 
a ] 3p | i (18) 
where P = p; — po > 0; and from Eq. (5) that 
2m > >3 
a =* (19) 
( p it 


If the independent variable in Eq. (17) is changed from ¢ to the volume ratio 


R3 
t= Bis = 2: ¥: (20) 
R - “s 
there results 
2 5 Id = 
i =~ “ees 2 i ee oe oe (21) 
ax oO 6 dx 6 


Equation (21) will be recognized as the differential equation for the hypergeometric 


function F(a, 8; y; x) where the parameters have the values 


—1 + 27(24n — 25)°° | na 
a= 19 = To 1 18; 
—1 — i(24n — 25)'” 
B = = as: (22) 
12 
l 
cae 


It is convenient to take the general solution of Eq. (21) in the form [6] 
a = AF(a, B:1/2;1 — x) + BO — 2)’ F(-—a + 1/3, —8 + 1/3;3/2;1 — 2 (23) 
where A and B are constants to be fixed by the initial conditions. If a = a, at ¢ = 0 or 


R = R, , then one has 
A = Gs (24) 


Similarly, if the initial velocity amplitude for the distortion is vp , then 


_— 8 ii E és 
im <  "—-., = dx dt_}’ 
= 8 | | 7 (25) 
2k, L3p 


° . fa bl . ‘ < 1/2 - — ° 
which fixes the constant B. The quantity (2P/3p)’° is a characteristic velocity and 
R, is a characteristic length for the system, and it is convenient to describe the initial 


from Eq. (23). 


velocity amplitude in terms of the length 
Volto 
ly = B53 > (26) 
(2F /3p) 
and Eq. (25) becomes 


B = 2l,/3. (27) 
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The limiting value of aas R >&, or x — 0, is readily determined from Eq. (23). One has 


‘R—4o) = = g'1T(9/ __% ee — — 
alk ) Qc x T(2, »| rd/2—a)- + 3) 17/6 +a | 
so that for large n 
-1/2_ tino / -1/6 l, <-7/6 
a. ~ ame" T(2 3)| a3 M645 | (28) 


where from Eq. (22) 
(24n — 25)'” 
12 


6= 


Figure 1 shows the variation of a/a) with R,/R for various values of n for the case in 
which the initial velocity amplitude is zero, l, = 0; Figure 2 shows the variation of a/a, 
for the case in which the initial velocity amplitude is different from zero. Of greater 
significance is the ratio of the distortion amplitude a to the mean bubble radius R; the 


— -- — —— —-, 
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EXPANDING BUBBLE WITHOUT SURFACE TENSION, "0 














Fic. 1. The ratio of the distortion amplitude a to its initial value ao is shown for an expanding vapor 

cavity as a function of Ro/R where Rp is the initial cavity radius and F its radius at later times. The 

distortion of the spherical interface is anY, where Y, is a spherical harmonic of order n. The initial 
velocity of the distortion is zero. 
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Fic. 2. The distortion amplitude is shown as a function of cavity radius for a non-zero initial distortion 
velocity. 
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behavior of (a/a,) (R./R) is shown as a function of R,/R in Fig. 3 for the case in which 
the initial velocity amplitude is zero, and in Fig. 4 for the case in which the initial velocity 
amplitude is different from zero. 
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EXPANDING BUBBLE WITHOUT SURFACE TENSION, 42:0 











Fic. 3. The ratio of the distortion amplitude a to the mean cavity radius F (in units of ao/Rpo) is shown 
as a function of R)/R for the case shown in Fig. 1. 
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Fic. 4. The ratio of the distortion amplitude a to the mean cavity radius FR (in units of ao/Ro) is shown 
as a function of Ro/R when the initial distortion ampiitude velocity is non-zero. 


(ii) Collapsing cavity; no surface tension 
Equation (21) is also applicable to this case. It is convenient, however, to write the 
solution in the form [6] 


a = Ax “Fla, a + 2/3; 1/2; 1 — 1/2) 
+ Bxr“(1 — 1/x)'°F(—B + 1/38, —B8 + 1;3/2;1 — 1/2), l<2<0. (29) 


Equation (18) is written in the form 


Hil _ 2p [z - | 
dt | 3p LR’ 
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where now p = ~» — p; > 0. The constants A and B are found as in the previous case in 
terms of the initial distortion amplitude a, and the initial velocity amplitude v, : 


A = Ao : 
2 Ro  _ _ 2 
B= 3 Vo (2p/: p)? 3 ly. 


The solution may, of course, be 


written in a variety of forms. In place of Eq. (29) one 
may write, for example, 


a = A’y*Fla, a + 2/3; 2a + 7/6; y) 


+ Bly *"'F(—a + 1/2, —a — 1/6; —2a + 5/6; y), (30) 


where now 


fa 
=ZzZ-—- Ss < < F 


A’ and B’ are linear combinations of a) and J, which will not be written explicitly here. 
From Eq. (30), one finds in the neighborhood of y = 0 that 


, 2+i - 2-16 
a= Aly 4 Big ie 
or 


a = const X R7'™*, (32) 


which is the singularity noted by Birkhoff. 


The variation of the distortion amplitude with mean bubble radius is shown in Fig. 


5 for n = 3 and in Fig. 6 for n = 6. The quantity of significance is the ratio of a to R; 


therefore, the variation of (a/ay) (R,/R) with R/R, is shown in Fig. 7 for n = 3 and in 
Fig. 8 forn = 6. 
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Fic. 5. The ratio of the distortion amplitude a to its initial value ao is shown for a collapsing vapor 
cavity as a function of R/Ro. 


The case shown is for a distortion described by a spherical harmonic of 
order n = 3. 
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Fic. 6. The distortion amplitude is shown as a function of cavity radius for the case n = 6. 
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Fic. 7. The ratio of the distortion amplitude a to the mean cavity radius R (in units of ao/Ro) is 


shown as a function of R/R» for n = 3. 


(iii) Expanding cavity with surface tension 
For a bubble expanding from rest, dR,/dt = 0, one has from Eq. (6), if surface tension 


is included, 
dR}? 2P R?] 20 R? ; 
@ | ~ 8p E " ze pR E i zi), (33) 


where P = p; — po > O. The radial acceleration is determined by the relation 
, aR 3/dR|/ P—2¢/R 
ki dt? + P |@ | i p 


If the stability equation (2) is written in terms of the independent variable 


Ro 
‘—s (34) 
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Fic. 8. The ratio of the distortion amplitude a to the mean cavity radius FR (in units of ao/Ro) is 
shown as a function of R/Ro for n = 6. 


it becomes 


2 = Y. 2 k 1 (3 ) ‘ da 
- ss “» i sees ~ ain wits om ao — - _ 2 
{2 ke + (i = 3} dz” {2 ~ 3x2 a dz 


—(n- nt [l — (n + 1)(n + 2)] + (1 — at) +h, = (, (35) 


where 


20 


k= Pp 


(36) 
so that k is the ratio of the initial value of the surface tension to the static pressure 
difference between the inside of the bubble and the liquid. Equation (35) has a neat 
solution for the special value of k = 2/3 in which case it reduces to the hypergeometric 
differential equation. This value of k is reasonable for vapor bubbles growing in super- 
heated water where it is effectively slightly smaller than unity [7]. A convenient form 
for the solution is 


a = AF(a, 8; 1/2;1 — 2) + BU — 2)'?F(a + 1/2, 8 + 1/2;3/2;1—2) (87) 
where now 


3, (9+ 16N)'” | 


a= ~ + 4 : (38a) 
« ( »aT\1/2 
8 es _3 = (9 + 16N he . (38b) 
4 4 
and 
wa? = Li (n® + 3n + 1). (39) 


If ao is the initial distortion amplitude and v, the initial distortion velocity amplitude, 
then one finds 
A= Ao , 
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and 
B = 2l,, 
where 
— Ro _ , 
. GPia”** 
The variation of a/a) with R,/R when k = 2/3 is shown in Fig. 9 for n = 2 and 3. 


In Fig. 10, the variation of (a/a,) (R,/R) with R,/R is shown for these same values of n. 
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Fic. 9. The distortion amplitude a relative to its initial value ao is shown for an expanding cavity as 
I aS - 


a function of Ro/R for the case in which the effect of surface tension is included. For n = 6 the curve 
with /)/ao = 1 is not shown since it lies quite close to the curve Ip/ay = 0. 
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Fic. 10. The ratio of the distortion amplitude a to the mean cavity radius (in units of ao/Ro) is shown 
as a function of Ro/R for the case in which the effect of surface tension is included. For n = 6 the curves 


lo/ao = 1 and Io/ao = O are very near each other. 
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When these curves are compared with Figs. 1 and 2, or Figs. 3 and 4, the stabilizing 
effect of surface tension is evident. It is also of interest to observe that, when surface 
tension is included, a/a,) changes sign as R increases. 

(iv) Cavity collapse under surface tension alone 

An additional case of special analytic simplicity occurs for p; = po so that the cavity 
collapses under the influence of surface tension alone. If the independent variable in 
the stability equation for a is changed from ¢ to the area ratio 


R? 
==) (4 
Re’ ® 


u 


the solution of the resulting differential equation is readily found to be 


a = u"[AF(a, 8; 1/2;1 — u) + BU — w)'?F(a + 1/2, 8 + 1/2;3/2;1-w], (4) 


where m has the value 


— (2 —_— 9ny'/s 
i ata : + en ee . (42) 


and F is the hypergeometric function with parameters a and 8 determined by the 


relations 


pa 2_©& — (n? + 3n + 4); 


a+B=2m+%. 


If the initial distortion amplitude is a, and the initial distortion velocity amplitude is 
vo , then 
A=Q, 
and 
B= -L, 
where 
Ry 


, (2¢ pR,)' . 


It is evident from Eq. (41) that 


a— const X R7'* as R— 0. 


3. Conclusion. For an expanding vapor cavity, an initially spherical shape is stable 
in the sense that the deformation amplitude a remains small compared with R if its 
initial value a is small compared with the initial cavity radius R, . The consistency 
and applicability of the linearized perturbation theory for the distortion amplitude is 
thus demonstrated. These conclusions from the linearized theory must be qualified for 
the case in which surface tension is negligible. As is shown graphically in Figs. 3 and 4, 
a/R as a function of R,/R has a maximum which increases slowly with n, the order of 
the spherical harmonic. It follows, therefore, when surface tension is unimportant, that 
needlelike irregularities in the spherical interface may grow to significant amplitudes. 
The present linearized theory is inadequate to follow the development of these high 
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order distortions of the interface. This instability for large n disappears when surface 
tension is of significance so that no such restriction need be imposed on the applicability 
of the linearized theory in this case. 

For a collapsing vapor cavity, on the other hand, the perturbation theory is valid 
provided the distortion amplitude is not followed to small cavity radii. If Ry is the 
initial radius of the spherical cavity, then the distortion amplitudes remain small so 
long as 1 > R/R, Z 0.2 where the lower limit is, of course, approximate. The linearized 
theory is thus applicable over an interesting and important range of cavity radius. As 
R — 0, the distortion amplitudes oscillate in sign with increasing frequency and increase 
in magnitude like R~'’*. This increase in distortion amplitude as R~'* is found with 
and without surface tension. It may be remarked that the linearized perturbation theory 
for the distortion amplitudes breaks down in a range of radii near that for which the 
present model of the vapor cavity becomes invalid. It is known [8] that the vapor pressure 
within a collapsing vapor cavity, such as is encountered in cavitating flow, begins to 
rise very rapidly as R/R, becomes smaller than approximately 0.1. 
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BY 
GARRETT BIRKHOFF 


Harvard University 


AND 
J. B. DIAZ 
Institute for Fluid Dynamics and Applied Mathematics, University of Maryland 


1. Flow problems. We shall be concerned with connected networks. These will be 
defined as finite connected graphs, on which the boundary is explicitly specified. 

As a graph’, such a network consists of a finite set N of nodes (or vertices), 
A, , ++: , A, , certain nodes being joined in pairs by a finite set L of oriented links (or 
branches) a, , --- , a, . Thus the graph is specified by an incidence matrix of n rows 
and r columns, |! ¢,; ||, where ¢,; is +1, —1, or 0 according to whether the node A, is 
the initial node, the final node, or not incident on the oriented link a; . It will be assumed 
that each link a; joins exactly two nodes, hence we may write a; = Aj <;)Ays,;) , Where 
A,;;) is the initial node of the link a; and A,,;, is the final node of the link a; . This 
implies that the incidence matrix has just two non-zero entries in each column (one 
being +1 and one —1). It will also be assumed that each node A, is incident on at 
least one link. This implies that each row of the incidence matrix has at least one non- 
zero entry. 

Further, a subset @N of N, called the boundary, is supposed to be specified. This 
subset 0N may or may not be empty. If dN is not empty then the elements of dN are 
called the terminals of the network. Finally, the network is supposed to be connected? 
in the usual sense that a graph is said to be connected. 

We shall consider first a special class of network problems, which we shall call ‘flow 
problems’. Whether they concern hydraulic networks or direct (electrical) current 
networks, flow problems involve two real valued functions: a potential function u(A;) 
defined on the nodes, and a current function 7(a;) defined on the oriented links. In hydrau- 
lie networks u(A,) is the pressure head; in direct current problems, it represents the 
voltage.** 

In network flow problems, leaks are neglected. One thus assumes, at each interior 
node A, in N — @N, Kirchhoff’s node law 


daa) =0, h=1,--- ,n; (1) 
j=1 


where, in view of the definition of the incidence matrix, the summation is effectively 


Received Nov. 1, 1954. 

*This work was sponsored in part by the U. 8S. Atomic Energy Commission. 

ISee W. H. Ingram and C. M. Cramlet [12], J. L. Synge [20], and the books by O. Veblen [7] and 
D. Kénig [9]. The basic ideas are due to G. Kirchhoff [1] and H. Poincaré [4, 5]. (Numbers in square 
brackets refer to the bibliography at the end of the paper). 

2Actually, this assumption plays a very small réle, but it simplifies the statement of various results. 

**Professor W. Prager has kindly drawn our attention to the occurrence of similar flow problems in 
the mathematical study of the distribution of traffic over a network of roads, see [28]. 
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taken only over the links incident on the node A, . For physical equilibrium, the currents 
must also satisfy certain equilibrium relations 


(2) 


i(a;) = c;(Au;), yg @ ],.+-< 5%, 


where Au; = u(A,;,;)) — u(A,,;)), with A,,;) the initial and A,,;, the final node, respec- 
tively, of the oriented link a, 

Physically, the conductivity functions c;(Au) are usually increasing and continuous. 
In our theorems below, we shall usually assume one or both of these conditions. For 


reference, we write” 


c;(Au) is a strictly increasing function of Au, (2a) 
c;(Au) is a continuous function of Au. (2b) 

Thus, in hydraulic problems, it is commonly assumed that 
c;(Au) = K;-sign (Au)-| Au |", where K,; > 0, a > 0: (2c) 


(For turbulent flow in pipes, a = 1.85 is commonly accepted.) In direct current problems, 
a linear relation 


i(a;) = c;Au; , (2d) 


is generally used. Since the general case will be considered in Sec. 4, we shall omit the 
physical condition c; > 0, which corresponds to (2c) with a = 1, and gives the classical 
case treated by Kelvin [2]. 

In summary, we will assume (1) at all “interior” nodes (i.e., the nodes of N — dN) 
and (2) on all links. At each node A, of the boundary dN, the total “influx” », must 


clearly satisfy 


>» 


y, = w €,;1(a;). (3) 
Comparing this last equation with (1), we get the necessary condition 
>» = 0, (3’) 


summed over ON. [This follows because 
> : ¢,,1(a;) = 0, 
A 1 


and (1) then implies 


>/ 


which is (3’).] 
To obtain ¢ 
terminal A, in ON; for example, one might assume 
I. The potential w(A,) is given, or 
II. The ‘total influx” v, at A, is given. 
Because of the obvious analogy with potential theory, we shall refer to a boundary 
value problem in which a condition of Type I is given at each terminal as a “Dirichlet 


1 “boundary value problem’, some condition must be given at each 


8The significance of (2a) and (2b) was first stressed by d’Auriac [14] and Duffin [17]. 
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problem’’. Similarly, if the total influx is specified at each terminal we shall speak of a 
‘‘Neumann problem’’. Problems involving both types have been treated in the literature‘. 
Still more generally, one can consider “mixed” conditions, of the type (notice that II’ 
includes II as a special case) 

II’. A functional relation », = F,(u) is given, where 


F’,(u) is a non-increasing, continuous function of u. (4) 


(In heat flow problems, this would correspond to a linear or non-linear “law of cooling’’.) 

2. Uniqueness theorem. It is not hard to prove a general uniqueness theorem, in- 
volving boundary value problems with boundary conditions of Types I, II, or II’, 
which is adequate for most physical flow problems. To formulate it, let w = u(A,) and 
i = i(a;); and wu’ and 7’ represent two different solutions of the same boundary value 
problem. Consider the expression 


D* = > (i, — i(Au, — Aut) 
L 


(5) 
= Zz: [e,(Au,) — c(Auj)][Au, — Au]. 
L 
In the linear case, clearly (5) simplifies to 
D* = > o,(Au, — Auli)’. (5’) 
L 


In any case, the following result is immediate: 
Lemma 1. If all the conductivity functions ¢,(Au) satisfy (2a), then D* 2 0. Strict 
inequality holds unless u — u’ is a constant. 
We now make a second evaluation of D*. By (3), 
Zz (c,(Au,) — c,(Aut)][Au, — Auf] = 7 {le,(Au,) — c,(Auz) }[ >> €,,(u(.A,) — u’(A,))]} 
J N 


L 


= > — vj)f{u(A,) — u’(A,)]; 


N 
where, for example, 


5, eee bis €,;1(a;) 
L 


is the “influx” corresponding to the potential u at the node A, of N. It follows from 
(5) and (1) that 
D* = >> (™ — vi)[u(A,) — u’(A,)]. (5+) 
aN 
For boundary value problems involving only conditions of Types I or II at the terminals, 
one has D* = 0. If conditions of Types I, II, or II’ occur, provided that (4) is assumed, 
clearly D* < 0. Comparing with Lemma 1, we get 
THEOREM 1. There is at most one solution to any boundary value problem defined 
by (1) and (2), with boundary conditions of Types I, II, or II’, provided that (2a) 
‘D’ Auriac [14] and Duffin [17] consider the Dirichlet and Neumann problems, plus a special “‘mixed”’ 
problem, where a Dirichlet condition is imposed at some boundary nodes and a Neumann condition at 
the remainder of the terminals. D’Auriac proves uniqueness and Duffin, existence and uniqueness 


theorems. 
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and (4) are assumed and that, in the Neumann problem, potential functions which differ 
only by a constant are considered to be identical. 

3. Dissipation function; variational principle for the Dirichlet problem. If 7 and u 
are any two functions defined on the oriented links and the nodes, respectively, of a 
network, we may define the dissipation function as the sum 


D = > i(a,) Au(a,). (6) 
L 


(The name “dissipation function” expresses the fact that, in the two physical problems 
mentioned in Sec. 1, the expression D represents the rate of energy dissipation.) We 
shall now derive an alternative formula for D, analogous to (5+) for D*. Since 


Au(a,) = u(Aiay) — WAsa) = D> exu(A,), 


one has 


dc (Au,)Au, = Zz. {c,(Au;,) > é,,u(A,)} = > y,u(A,); 
L N 


L N 


and thus 
D = >> »u(A,). (6+) 
ON 
The expression D, according to (6+), represents the rate of energy influx. 

In the linear case, the dissipation function reduces to pe c,(Au,)’, and it is classical® 
that this is minimized by the solution of the network problem over the class of potentials 
assuming the given terminal potentials. We shall now derive an analogous variational 
principle for the non-linear case. However, this will not, in general, involve the dissipa- 
tion function. 

To formulate the new variational principle, we suppose u is given on dN, but is 
unknown on N — ON. For any assumed values of u on N — ON, we can then satisfy 
(2) automatically by defining i(a,) = c,(Au,) on each link a, . It remains to satisfy (1), 
and for this we shall find a variational formulation. Namely, define the functions C, by 


pAu 


C (Au) = | c,(x) dx, 
70 


so that the derivatives 
dC, 


Ci(Au) = ~-= = c,(Au); 
z( Au d(Au) c,.( Au) 


for simplicity, we shall assume condition (2b). (Duffin [17, pp. 965-967] uses this same 
device of auxiliary functions for what we call the Neumann problem.) 

THEOREM 2. For given u(A,) on ON (i.e. for the Dirichlet problem), assuming (2b), 
the first variation of 


Vu) = >> C,(Au,) (7) 
L 


is zero at each (“interior”) node of N — @N if and only if Kirchhoff’s node law (1) 
holds at each interior node. 


‘See W. Thomson [2]; J. C. Maxwell [3, vol. I, pp. 403-408]. 
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Proof: (By the first variation is of course meant the following limit: 
| 
s oa 
6V(u) = — Viu+ edu)!| , 
de le=0 


where 6u is any potential function defined on N but which vanishes on dN.) By direct 
computation, writing a, = A; «Ayu , one has 


D> Ci(Au,)[6u(A ia) — bu(Asa)] 
L 

D> {Ci(Auy) px €,,6u(A,)} 

L N 

Dd {6u(A,) Do eni(a,)}. 

N L 


For A, in dN, the number 6u(A,) is zero, while for A, in N — ON, the éu(A,) are arbi- 
trary. The conclusion of Lemma 2 is now evident. 

Coro.uuary. If (2a), (2b) hold, then Kirchhoff’s node law (1) holds if and only if 
V'(u), considered as a function of the arbitrary values of the potential at interior nodes, 


5V 


I] 


Il 


II 


has an absolute minimum. 

For if, regardless of (2a), the function V(u) has even a local minimum, then 6V = 0, 
and hence Kirchhoff’s node law holds. While, on the other hand, if (2a) and (2b) hold, 
then V(x) is a convex function of u, since it is a sum of convex functions either of the 
individual u, = u(A,), or of pairs of these variables, as may be readily seen from (7). 
We leave the detailed verification of this to the reader. Hence, if Kirchhoff’s law holds 
for some u, the convex function V(u) must have an absolute minimum for this particular wu. 

Remark. In the case (2c) of an exponential conductivity law, with the same exponent 
a for all links in the network, the dissipation function D is proportional! to the function V, 
and therefore D can be used in place of V in the results of this section. 

4. Existence theory for the Dirichlet problem. We shall now derive an existence 
theorem for the Dirichlet problem which is adequate for most physical applications. 
In order to avoid giving the impression that it is the “‘best possible”, we shall preface 
it by giving a much stronger result for the linear case. 

In the linear case (2d), a given trial potential function, when used to construct 7 by 


means of i(a,) = ¢,[u(A;a)) — u(A;a))], for each link a, = A; «)Aysau) , Will satisfy 
Kirchhoff’s node law (1); i.e., (see Sec. 3) will solve the Dirichlet problem, if and only if 
De ente(U(A cc) — WAray)] = 0, (8a) 
k 
for every A, in N — @N. This gives a system of s linear equations in the s unknowns 
u(A,) = u, , which may be more compactly written thus 
> e,;u(A,) = b,, h=1,--- ,8, (8) 
j=1 


where the numbers b, are known. The matrix of coefficients || ¢,; || of (8), which is sym- 
metric (as follows readily from (8a) and the definition of the incidence matrix || €. ||) 
will be called the conductivity matrix of the network. It is well known® that for any 
system like (8), existence and uniqueness are equivalent to each other, and also to the 


6G. Birkhoff and S. MacLane [25, chap. X). 
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condition that the determinant of the matrix of coefficients be different from zero. 
This gives the following result’. 

THEOREM. If (2d) holds, the Dirichlet problem is solvable, for a given network N 
(with at least one interior node and at least one boundary node) for arbitrary values, 
if and only if det || ¢,; || + 0. This condition is also necessary and sufficient for uniqueness. 

We now see how special, in the linear case, is the condition (2a) requiring all con- 
ductivities to be positive. If this condition holds, then all the diagonal elements ¢,, 
are positive, while 

Ci, = >> Ci for h =i, ---,28, 


1A 


with strict inequality holding if and only if the node A, is linked directly to a boundary 
node (which will certainly occur for at least one node, since neither N — ON nor ON is 
empty, and the network is connected). It follows” that, if, in addition the matrix || c,, 


is not reducible to the form 
P U 
0 Q 


by the same permutation of the order of the rows and columns, where the matrices 
P, U, Q, 0 are all square matrices, and 0 consists only of zeros, then det || c,, || # 0. 
However, it is easy to see, again using the theorem mentioned in footnote 8, that if all 
the conductivities are positive, and the conductivity matrix of a connected network has 
the ‘‘exceptional” form just mentioned, then its determinant is still not zero. For if 
'! c,, || is of this exceptional form then its determinant is the product of the determinants 
of P and Q, each of which is again symmetric and “dominantly diagonal’, and may be 
further reduced, in the same way that the original conductivity matrix was reduced, 
in case either of them is exceptional. (Notice that, in view of the symmetry of |! cy, ||, 
it follows that the submatrix U must consist only of zeros.) Continuing this reduction 
as far as possible until only non-exceptional symmetric matrices occur (only a finite 
number of steps are possible) one finds that the det || c,, || is the product of a finite 
number of determinants, each corresponding to a “dominantly diagonal’ matrix which 
is not “‘exceptional’’, and that all elements not appearing in this product are zero. Since 
the given network is connected, at least one node in each subnetwork associated with 
these submatrices must be linked directly to a boundary node of the given network. 
Hence, by the theorem mentioned in footnote 8, the determinant of each subnetwork is 
not zero, and thus det || ¢,, || is not zero either. It is clear that this class of non-singular, 
dominantly diagonal symmetric matrices is but a very small subclass of the class of all 
non-singular symmetric matrices. 

The existence theorem which we shall now prove for the (possibly) non-linear case 
corresponds, however, to the theorem (in the linear case) obtained from Theorem 2 
upon making the superfluous additional assumption that the conductivity matrix 

c,; || is “dominantly diagonal” in the sense just described above. 


7See C. Saltzer (27, p. 122], J. L. Synge (20, p. 127). 
8See Theorem III of Olga Taussky [18, p. 673]. For an application to electrical networks, see M. 
Parodi [13]. 
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By Theorem 2, any local minimum of V(x) will provide a solution, in the non-linear 
or linear case. However, if every C,(Au) — +o as | Au | — +o, then V(u) will be 
bounded below everywhere; and be arbitrarily large’ outside any sufficiently large 
bounded “cube” in (wu, , --+ , u,) space. Hence V(u) will have an absolute minimum 
inside some such cube, by a theorem of Weierstrass on continuous functions. We conclude 

THEOREM 3. If (2b) holds, and if, for all k, 


r ¢.(x) dx = [ , c,(x) dx =+o, (9) 


in the sense of improper Riemann integration, then the Dirichlet problem has a solution 


for arbitrary boundary values. 
Coro.uary. If (2a) and (2b) both hold, then (9) may be replaced by the conditions 


c.(x) > 0, for some x> 0, (9a) 


and 
c,(x) < 0, for some s < ©. (9b) 


5. Neumann problem. The Neumann problem is dual to the Dirichlet problem, in 
the sense that the réles of u and 7 are interchanged. To make the duality more marked, 
we note that, since any continuous, strictly increasing function y = c¢,(x) has a (unique) 
continuous, strictly increasing inverse function x = r,(y), conditions (2a), (2b) are 
self-dual. Accordingly, we shall replace (2) in Sec. 1 by 


Au, = r;,{i(a,)], (10) 


and refer to the r, as resistance functions. The condition that there exists a single-valued 
potential u(A,), such that Au, = u(A;) — u(A;) whenever a, = A,A;j , is evidently 
Kirchhoff’s circuit law 
> r,[i(a,)] = 0, (11) 
; 
for any sequence I of oriented links forming a closed cycle (or circuit). 

For a given influx v on AN, satisfying the consistency conditions (3’), the most 
general current function 7 which satisfies Kirchhoff’s circuit law (11) is obtained by 
‘adding’, onto some fixed current function satisfying the same conditions, “cyclic” 
currents 8, , , 8, around closed cycles T, , --- , I, forming a basis for the closed 
cycles of the network. This fact is easily seen in the case of a planar network (graph), 
when the basic cycles may be taken as the (oriented) boundaries of the polygons into 
which the network subdivides the plane, and one has r + 1 = n + ¢. The general case 
is also elassic’® (i.e., 4 = r — n + 1 for a connected graph). 

Thus, once an initial current distribution satisfying (3’) and (11) has been found, 
each 6 = (8, , --- , 8.) determines a unique current distribution on the set of links L, 
satisfying (3’) and (11), while (10) may be taken as defining Au. (We treat (10) as a 
substitute for (2), recalling that (2) and (10) are equivalent if (2a) and (2b) hold.) 

We now define R,(z) = fo r.(y) dy, for each link a, , so that the derivative Ri(z) = 
dR,/di = r,(z), and assume for simplicity that the r,(y) are continuous. 


°This follows from a modification of an argument of Duffin [17, p. 965] which uses in an essential 


way the fact that the network is connected. 
Poincaré [4], Veblen [7, p. 9], Ingram and Cramlet [12, p. 137], and Synge (20, p. 123). 
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THEOREM 2’. For given consistent values [see (3’)] of », on ON, the first variation 
of the function 


W(s) = > R,{i(a,)] (12) 
L 


vanishes identically if and only if the r,[z(a,)] satisfy Kirchhoff’s circuit law (11). 
Proof: By direct computation, 


bW = >> Rifi(a, éi(a,) = >> 68; >> r, fia], (13) 

E B T ; 
where the last sum is taken over a basis B of the closed cycles of the network, which 
consists of T, , --- , ', . This last sum is zero for arbitrary 6 if and only if the individual 


sum taken over each I; is zero, which is equivalent to (11). 

Coro.uaRry. If (2a), (2b) hold, then Kirchhoff’s circuit law (11) holds if and only if 
W(8) has an absolute minimum. 

The proof is similar to that of the corollary to Lemma 2, and may be omitted. 

Remark. In the case (2c) of an exponential resistance law, with the same exponent 
a for all links in the network, it follows (see Sec. 3) that the dissipation function D is 
proportional to the function W, and therefore D can be used in place of W in the above 
results. This is known in the linear case’’. 

THEOREM 3’, The Neumann problem has a solution for any set of compatible bound- 
ary influxes [see (3’)], provided that (2b) holds and that 


[ r.(y) dy = | rdy) dy =+o. (14) 

The proof is identical with that of Theorem 3, and the analogue of the corollary to 
Theorem 3 also follows similarly. 

6. Mixed boundary value problem; relaxation methods; existence theorem. Without 
striving for maximum generality, we shall prove an existence theorem which is adequate 
for most applications. The method of proof to be employed is constructive, in that in 
many concrete instances the performance of the “relaxation steps” used in the proof 
can actually be used in order to construct numerically a solution to a network boundary 
value problem. 

Let us suppose that our boundary conditions are of Types I and II’. Since Kirch- 
hoff’s node law (1) really corresponds to a condition of Type II’, with F,(u) = 0, we 
can reformulate our problem as that of satisfying a condition of Type II’ at all those 
nodes A, , --* , A, where the potential u(A,) zis not prescribed; we shall denote this 
set (supposed to be not empty) of nodes by M. Further, we shall suppose that the set of 
nodes at which the potential zs prescribed, which is N — M, contains at least one node. 
It is for this class of boundary value problems that an existence theorem will be proved, 
under the assumption that (2a) and (2b) hold and that 
lim ¢,(2) =+o. (15) 


z—++@ 


lim ¢,(x) =— ©; 
z—7—@ 
Thus we shall exclude the case of ‘“‘saturation currents’’. 

(The requirement that the set VN — M be non-empty, seemingly—but only seemingly 
—excludes the Neumann problem from consideration. Because, granting, for the purposes 





11W. Thomson (Lord Kelvin) [2], and J. C. Maxwell [3]. 
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of the present discussion, that an existence theorem has been proved for the above 
mentioned class of mixed problems, then an existence theorem for the Neumann problem 
readily follows from it. This can be seen by merely assigning arbitrarily the value of 
the potential at a fixed node of the network, as an additional boundary condition, 
besides the given Neumann conditions. By this obvious artifice, any Neumann problem 
can be turned into a mixed problem of the class described above, and hence has a solution 
for each arbitrarily assigned value of the potential at the chosen fixed node, i.e. a “‘one- 
parameter” family of solutions. In view of this, the Neumann problem need not be 
mentioned in the following discussion.) 

Just as in Sec. 3, we can satisfy (2) by fiat for any choice of u, = u(Ay), «++, Um = 
u(A,,), merely by defining z(a;) = c;(Au;) for each link a; . We can then compute 
v»% = >-1 &;%(a;), for each node A, in M, and define the discrepancy (or residual) function 


6, =», — Flu), h=1,-+-+,m. (16) 
An existence theorem clearly asserts that 6(u) = 0 for some u = (u, , -++ , Um). 
Lemma 4. The function 6 = TJ (u) is one-to-one and continuous. 


Proof: The continuity of T follows from the fact that the functions c, are continuous 
by (2b), and the functions F’, are also continuous, by (4). It remains to show that distinct 
u determine distinct 6. This follows readily from Theorem 1, but we shall go over the 
proof, to emphasize the réle of the requirement that the set N — M, where the potential 
values are assigned, is not empty. To this end, consider, as in (5) and (5+) that 


D* = p (a, — ti) (Au, -— Aut) > 0, 
L 


by (2a). On the other hand 


D* = DO {(e — tf) Do enlu(A,) — w(A,)]} 
dD (vs — v4)[u( Ar) — w’(A,)], 


p> (v, — vf)[u(A,) — w’(A,)] 


and if T(u) = 6 = 6’ = T(u’), then 
D* = > (Film) — F’(w)][u(A,) — w(A,)] 
M 


IIA 


0, 


by (4). Hence D* = 0, and by Lemma 1, it follows that u — wu’ is a constant. But this 
constant difference must be zero, since it is zero for each node in N — M, which is not 
empty. 

Now, still assuming (2a), (2b) and (15), we pass on to a relaxation method. We 
shall consider the residuals” 6, of a variable trial function u(A,), which are defined by 
(16). We first prove four lemmas involving the “order” relation. 

Lema 5. As u(A,) is increased, all other values of u being held fixed, 6, increases, 
all “adjacent’”’ 6, decrease, and all other 6; remain constant. 

Proof: For if A, is “adjacent” to A, , that is, there is either a link A,A, or a link 
A,A, in the network, then an increase in u(A,) increases [by (2a)] either 7(A,A,) or 


122We shall conform to the terminology of R. V. Southwell [11], where possible. 
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—i(A,A,), as the case may be; hence it increases v, , decreases vy, , and leaves unchanged 
vy; when A; is not adjacent to A, . Also, by (4), an increase in u, = u(A,) either decreases 
or leaves unchanged F’,(u,), and leaves unchanged all remaining F;(u;), where A; # A, . 

Lemma 6. Consider a node A, , and suppose that the values of wu at all adjacent 
nodes A, are increased, while the values of u at A, , and at all nodes not adjacent to 
A, are held fixed. Then 6, decreases, while 6, , where A, is adjacent to A, , increases. 

The proof follows along similar lines to that of Lemma 5. 

Now, consider 6, [see (16)] as a function of the single real variable w(A,), all the 
other u(A,) being kept constant. From (2b), (15), and (4) it follows that 6, is a strictly 
increasing continuous function of u(A,), which varies continuously from — © to +© 
as u(A,) does the same. Hence there is exactly one choice of u(A,) which will make 
6,[u(A,)] = 0 (“liquidate the residual’ at A,), provided that all the other w(A,) are 
kept constant. We define (exact) point relaxation at each node A, to consist of replacing 
the value u(A,) by this particular value which makes 6, vanish at A, , all the other 


’ 


values u(A,), for A, # A, , being kept constant. 

Lemma 7. Relaxation at a given node is isotone on trial solutions of the same problem 
(i.e. it preserves order). 

Proof: The proof is by contradiction. Suppose that u and w’ are such that u(A,) 2 
u'(A,) for all A, in N, and let v, = v(A,), and vj = v’(A,) denote, respectively, the 
functional values obtained from u and wu’ by point relaxation at the node A, . Suppose, 
contrary to what we wish to prove, that v, < v; . Now, starting with the “relaxed” 
function uZ (i.e. with the function whose value at each node A; # A, is u’(A;), while 
at A, its value is vs) one can proceed in two steps to the “relaxed” function up, , and 
obtain a contradiction, as follows. First, replace the values of wu’ at all the nodes different 
from A, by the corresponding values of u at these nodes, leaving the functional value 
unchanged at A, itself, and denote the resulting ‘“‘hybrid”’ function u’*. By Lemma 6, 
and the definition of point relaxation, one has that 


0 = 6,(uz) = 6,(u’*). (17) 


Secondly replace the value of u’* at A, , which is v; , by v, , and leave the values of u’* 
at all nodes different from A, unchanged. The resulting function is precisely the “‘relaxed”’ 
function uz . Since, by assumption, v, < v; , it follows from Lemma 5 that 


5,(u’*) > 6,(up). (18) 


But a comparison of inequalities (17) and (18) then shows that 6,(ug) < 0, contradicting 
the fact that, since v, was obtained by point relaxation of u at A, , the number 6,(upg) 
must be zero. This completes the proof of Lemma 7. 

In the proof of the following lemma and the theorem to follow we shall make use of 
two more additional assumptions, one concerning the conductivity functions c, and 
the other concerning the functions F,, of (4). For convenience we write them as follows: 


For every k, one has ¢,(0) = 0, (19) 


For every function F’, which does not vanish identically, there is a number 
x, such that F(x,) S 0 and a number 2, such that F(z.) 2 0. (20) 


In view of (4), it follows from (20) that whenever F, is not identically zero then it is 
= 0 for all sufficiently negative x and that it is < 0 for all sufficiently positive x. As for 
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(19), it certainly holds in the important special cases (2c) and (2d), and it means intui- 
tively, that “if the potential is constant then there is no flow of current”’. 

LemMMA 8. Suppose that (19) and (20) hold, in addition to (2a), (2b) and (15). 
Let uw be an arbitrary trial function (i.e. having the prescribed values on N — M),. 
Then there exist two other trial functions v) and w, such that 


vo(A,) < Uo( A,) < W,(A,), 


and 
5,(¥) < O SF &,(wo), h=1,-++-,™m. 


Proof: It will suffice to show how to construct the trial function v» such that both 


v,(A,) S U(A,) (21) 

and ; 
6,(v%) S 0, (22 
forh = 1, --- , m, since the construction of wy, is entirely analogous. The function v9 


will be defined in the following manner: 
\c, for A, in M (23) 


v,( A,;) = 
lnfA). t <AmN~WM, 


where C is a constant, which is to be chosen sufficiently negative so that the require- 
ments (21) and (22) asked of v) are met. First of all, if C S miny uw, then (21) clearly 
holds. As for (22), notice that if : in M is not linked to any node of N — M, and F, = 0, 
then [by (19)] it follows from (23), for any choice of C in (23), that 5,(vo) = (vo) = 0, 
and (22) holds; however, if F, 2 0, then still »,(v)) = 0, so that [by (20)] by choosing 
C sufficiently negative it will be true that 6,(v.) = —F,(vo) < 0, and (22) will again 
hold. It remains to consider the case when A, in M is linked to at least one node in 
N — M. From (3), in view of (19), it follows that, for any choice of C in (23), only the 
links joining A, to a node of N — M contribute essentially to the sum in »,(v9), and 
from (2a), (2b), (4) it is seen that if C = »)(A,) is sufficiently negative, then 6,(v)) = 
v,(Uo) — F (vo) will be S 0, fulfilling (22). Thus, all in all, in order that v, defined by 
(23) fulfill the requirements (21), (22), one sees that C must satisfy a finite number of 
conditions, all of which may be made to hold simultaneously, if only C is chosen suffi- 
ciently negative. This completes the proof of Lemma 8. 

We are now ready to prove our main result’* 

THEOREM 4. Suppose (2a), (2b), (15), (19), (20) hold. Let uw» be any initial trial 
solution of a mixed network flow problem, and suppose u, , U. , Us , --* are obtained by 
successively “‘point-relaxing”’ the residuals of the initial “ag solution up at an infinite 
sequence of nodes of M, in such a way that each node A, in M occurs infinitely often 
in the sequence of nodes. Then the sequence of trial functions u, , U2. , Us , +++ converges 
to the solution z of the given problem, the uniqueness of which has already been estab- 
lished in Theorem 1. 

Proof: First, by Lemma 8, there exist trial functions vy and wy such that both 


Vol 2 1,) = Uo( 2 4,) s Wo(A,), 


and 
5,(v) SO S 6,(wo), h=1,-+++,m. 


8This generalizes directly a result of J. B. Diaz and R. C. Roberts [22 
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Let v, , u; , W; denote the functions obtained from v , Up , Wo , respectively, by point 
relaxation at the first node of the preassigned sequence of nodes. By Lemma 7, the 
initial sandwich order is preserved, i.e. 


v,(A,) S u,(A,) S w,(A,), Rim 1, +++ mM. 


As a matter of fact, since 


6,(vo) S O SF 5,(w»), hw=i,-+* .m, 


it actually follows that (see Lemma 5) 


vo(A,) S v,(A,) S u,(A,) S w,(A,) S w(A,), h=1,---,m 


and that 

6,(,) S$ 0 S 6,(w,), h=1,-+> ,m. 
Similar inequalities hold for any positive integer n, if we denote by »v, , u, , W, , respectively, 
the functions arising from v , Uo , Wo , respectively, after successive point relaxation at 


the first n nodes of the preassigned sequence of nodes. Namely, we have 


Vo( A,) < v,(A,) S.s¢ = v,(A,) 4 u,(A,) = w,,(A,) ae (24) 


<= w,(A,) S w,(A,), 
and 
5,(v,) £ 0 S 5,(w,), homi1,-+++,m. 
Since, for each A, in M, the sequence of numbers v(A,), v,(A,), -*- , Un(As), «°° 18 
non-decreasing and bounded above [e.g., by w(A,)] it follows that the following limit 


exists 


v(A,) = lim2z,(A,), fim ll e<- om. (25) 
For A, on N — M we have that v(A,) equals u.(A,), which is exactly the value each 


function v, has at A, . Thus, to show that the function v is indeed a solution of the mixed 
problem, it only remains to show that 


6,(v) = 0, fim dl o+- om. (26) 
To do this, consider A, in M. Since A, occurs infinitely often in the preassigned sequence 
of nodes employed in point relaxation, it follows that there is an infinite sequence of 
positive integers n, < nm. < mn; «++ such that 
5,(v,,) = 0, k= 1,2,3,--- 
But then from (25) and the continuity of 6 (see Lemma 4), Eq. (26) follows. 


By proceeding in a similar manner with the non-increasing, bounded below sequence 
of numbers w,(A,), w,(A,), «++ , Wa(Aa), «++ one obtains that the function w defined by 


w(A,) = lim w,(A,), (27) 
for A, in N, is also a solution of the mixed problem. The uniqueness Theorem 1 then 


shows that v = w, the solution of the mixed problem and finally (24), (26), and (27) 
then show that u, also converges to the solution of the mixed problem. 
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—NOTES— 


UNSTEADY VISCOUS FLOW IN THE VICINITY OF A STAGNATION POINT* 
By NICHOLAS ROTT (Cornell University) 


Consider a steady two-dimensional “‘stagnation point”’ flow of a viscous incompressible 
fluid in the upper z, y-plane; let the flow be directed towards and limited by a plate in 
the plane y = 0, with the stagnation point at x = y = 0. The corresponding flow pattern 
is well known as an exact solution of the Navier-Stokes equations. 

Now, in addition, let the plate perform a harmonic motion in its own plane, i.e., 
in the x direction, while the flow at y — remains steady. It seems to have remained 
unnoticed that even in this case, the exact Navier-Stokes equations yield a soluble 
problem of the boundary-layer type. Using for the velocity components u and »v in the 
z- and y-directions respectively, 


u = axf'(n) + be'°'g(n), (1) 
v = —(av)'’’f(n), (2) 
and for the pressure, p, 
p= —* ax — pvaF(n) + po , (3) 
where 
7 = (2). (4) 
Vv 
(v = kinematic viscosity), and introducing these expressions in the Navier-Stokes 
equations, the following set of equations is easily obtained: 
f? —ff"’ =14+f", (5) 
kg + of’ — fg’ = 9", (6) 
ff’ = F’ — f”, (7) 
where k = w/a is a “reduced” frequency. Equation (6) and Eq. (7) result from the 


equation of motion in the z-direction, by putting the terms proportional to x and in- 
dependent of z respectively equal to zero. It is seen that Eq. (5) for the steady part is 
independent of the superimposed “g-flow.’”’ With the boundary conditions, f(0) = 
f’(0) = 0 and f’(~) = 1, Eq. (5) has the well-known Hiemenz solution. The viscous 
pressure term F can also be computed independently of Eq. (6). 

Since the function f(n) is known, Eq. (6) can be solved for g with the boundary 
conditions g(0) = 1, g(~) = 0. Consider first the steady motion of the plate with con- 
stant velocity b, i.e., # = k = 0. The corresponding solution gp , say, fulfills the equation, 


go + f9o — f'go = 9. (8) 

*Received Nov. 12, 1954. This research was supported by the United States Air Force, through 
the Office of Scientific Research of the Air Research and Development Command. The author wishes 
to acknowledge helpful discussions with Professor C. C. Lin during the latter’s stay at Cornell University. 
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The exact solution is 


es. ” 
i = f’(0) = Silf (9) 


or, the velocity profile is proportional to the shear distribution of the Hiemenz flow 
(see Fig. 1). Solution (9) obviously fulfills the boundary conditions since f”’(“) = 0, 


: — 


Jo 




















E ~ 3 


Fic. 1. Steady and quasi-steady velocity profiles. 


0 





and also the differential equation (8), because upon differentiating Eq. (5), the result is 


_— +4 f a a = 0. 


With f’’(0) =—1, the shearing stress at the wall, 7» , is proportional to 
f’""(0) 1 
(0) = * laa ~~ "ait 9 (10 
om a 


so that the resulting +» from both the f- and the g-flow becomes 


tw = Ha)"\ax5"0) a ath (11) 


Note that for « = b/a, the velocity outside the boundary layer is zero relative to the 
wall; the shearing stress at the wall, however, is zero at x = .658b/a. Evidently ry = 0 
does not necessarily mean separation for moving walls; in a system where the wall is at 
rest, the flow is not steady. The resulting boundary layer profiles are sketched in Fig. 
2. Such development may be expected at the stagnation point of a Flettner rotor. 


~~ 





-1 0 — 1 ax/b 


Fic. 2. Boundary layer development in the vicinity of a stagnation point with the wall moving at 
constant velocity b. 
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In the oscillating case, the two limiting cases, k < 1 and k > 1 will be considered. 
For k < 1, put 
g = go + tkg, + (tk)’g2 + ---. (12) 
Equating powers of zk, the solution for go is given by Eq. (9), obtained above; g, fulfills 


the equation 


gn’ + foi — f'n = G = B1lf”’ (13) 
with the boundary conditions g,(0) = g,(~) = 0. One particular solution of the in- 
homogeneous equation (13) can be immediately given: g, = .811f. This function, how- 


ever, does not fulfill the conditions at infinity, as f(“) — 7. The complete solution of 
Eq. (13) is obtained by adding the general homogeneous solution: 


f() ° [ E : 
n= sr tef” | an af” ) 
Ji f (0) + ¢;; i f 3 dn + Cof ’ (14 
where 
n 
E =e (~f f dn). (15) 
From g,(0) = 0, it follows that c. = 0. The constant c, has to be adjusted so 
that g:() = 0. Its value can be found with the help of the identity 
, ane 3 od inset 
f=f" i rads dn J, r dn (16) 


which may be proved by identifying the proper inhomogeneous solution of Eq. (13) 
with f, or directly by using the properties of the function f as a solution of Eq. (5). It 
can be shown also that for 7 ~~, f’”/E — 0, and that the integrand of the inner integral 
in Eq. (16) is integrable between the limits 0 and ©. Therefore, for large n, the inner 
integral in Eq. (16) will vary only slightly and may be replaced asymptotically by the 
constant value obtained after taking the limits between 0 and ~: 


nod fir? ) an E 
j sa tf E dy “mg yir2 dn. ( 16a) 
‘ JO E f- 0 J 
With this asymptotic expression, c, can be determined, and the final result is given by 
ee ae 
a <a —— “— t =s dn. 
A=F'Q- WO) EU J, 77 (17) 


The curve g;(7) is also plotted in Fig. 1; its extremum is —.152. The shearing stress at 
the wall is proportional to 


l 


TOF | - dn = —.496 (18) 


gi(0) = 


so that 7» due to the g-flow is, to the first order in 7k, 
tw = —plav)'’*b(.811 + .496ik)e***. (19) 
For high frequencies, k >> 1, the WBK-method is appropriate. In Eq. (6), put 


g = exp ( 8 dn) (20) 
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then, s has to fulfill the equation 


s’ +s’ + fs — f’ = tk. (21) 
Now set 


5. = (ik) '/*8 +s, + (ik)~'”*s, + (ik) ~'ss + +e (22) 


which is also equal to 


k 1/2 1 1/2 1 1 1/2 
en -(5) iadiealal (si) ss +i (5) me 
{lk 1/2 1 1/2 1 1 1 1/2 . 
- d(§)"s - (Ge) +p - EGR) hee 


Introduction of Eq. (22) into Eq. (21) gives upon equating powers of (ik)'”*: 
s=1, =f, & = af + if, (24) 
8s = —bff’ — if”,-:-. 


The boundary condition g(0) = 1 is always fulfilled by the choice of limits in the integral 
in Eq. (20), and g(~) = 0 is assured by the selection of the proper sign if the double- 
valued quantity (ik)'’’, namely, the one indicated in the decomposition Eq. (23). The 
resulting profile, using 8 and s, only, is 


y= eo[-()"a] -o[-)"s]oo(-fs4) os 


The first two factors represent the “Stokes-solution” [1], which would be obtained by 

complete disregard of the f-flow. The last factor has the value of about .92 for » = 1 

and .59 for » = 2. If, for large values of k, the second factor in Eq. (25) already assures 

a strong decay with 7, the last factor remains very close to 1 in the whole domain where 
significant values of g may be found. 
The shearing stress at the wall is 

tw = p(av)'’*bs(O)e'*'. (26) 


Now the value of s(0) differs from the Stokes value only if terms up to s; are included. 
From Eqs. (24), Eq. (23) gives 


«- -()"-{("- gre] * 


Therefore, the ratio of the out-of-phase shearing stress, tw, , to the in-phase component, 


Tw, ,18 


“$i = 1 — 654k”, (28) 
Tw, 

whereas for low values of k, from Eq. (19), the ratio is 

“Ft = 612k. (29) 

TW, 


In Fig. 3, the limiting cases for large and small k, Eqs. (28) and (29), are plotted, 
together with an estimated curve which joins the two plots smoothly. It is seen that for 








448 NOTES [Vol. XIII, No. 4 





ia h 4 ee 
4 pe on 
_ oe PIKE 


VA 





ZA 
ws 

















0 1 2 k 3 


Fic. 3. The ratio of the out-of-phase and the in-phase components of the shearing stress at the wall 
as function of the reduced frequency k, for oscillations in z-direction (g) and z-direction (h). 


a reduced frequency k larger than 2, the Stokes’ value rw,/tw, = 1 gives errors less 
than 20%. 

For the case of complex w, or exponential acceleration, only one special case will 
be considered, for which the solution may be obtained without any computational 
effort. Put k = — i, i.e., let the time-dependence of the plate velocity be given by the 
factor e*’. The corresponding special equation for g is 

gq’ + fg’ —(f’' + Dg = 0 (6a) 
which has the solution 
(30) 
in view of Eq. (5). The complete solution is 


u = axf'(n) + be"[1 — f'(n)] 


and the wall shearing stress is 


tw = play)’ f’’(0)[ax — be*‘]. (31) 
Here it is seen that, in contrast to the case k = 0, Eq. (11), the zero of the shearing 


stress and the zero of the relative velocity between plate and fluid occur at the same 
value of 2[=(b/a)e*']. 

Now let the plate move in the z-direction perpendicular to the xy-plane, i.e., the 
plane of flow. If the plate motion is uniform, the case under consideration is identical 
with the problem of the stagnation point in yawed flow, which has been solved by 
Prandtl [2] and Sears [3]. For unsteady motion of the plate in z-direction, Wuest [4] has 
already pointed out that the problem is soluble, and has given some numerical examples. 
Both cases become exact solutions of the Navier-Stokes equations, if the basic steady 
flow is the stagnation-point flow in the half-plane, as is presently assumed. If w, the 
velocity component of the flow in the z-direction, is taken in the form 

w = ce'”'h(n), (32) 
then Eqs. (5) to (7) remain unaffected, and the Navier-Stokes equation in z-direction 
yields 

ikh — fh’ = h’’ (33) 


with the boundary conditions h(0) = 1, h(~) = 0. For w = k = 0, Prandtl’s solution 


is obtained: 
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[ E dn 
ho = 22 (34) 


[ E dy 


[see Eq. (15)]. A series development for small k analogous to Eq. (12) leads to the equa- 





tions 


hi’ + fh, = 0, hi’ + fhi = ho , hi’ + fh; = h,, °°. (35) 
The solution of the second equation, adjusted to the boundary conditions h,(0) = h,(“) = 
(Q), is: 


—s [ (1 — ho)ho dn + (1 — ho) [ % dn. (36) 


7 
/0 ho to 


The profiles h, and h, are plotted in Fig. 1. The shearing stress at the wall in the z- 
direction becomes, to a first approximation in k: 
tw = —plav)'*c(.571 + .685ik)e'*'. (37) 


The investigation of the flow-behavior for k > 1 follows completely the method 
already employed for the g-flow. Putting 


h = exp (/ r dn). (38) 


in Eq. (33) yields 


re +r + fr = tk (39) 
and the series development analogous to Eq. (22) gives 
ro = 1, r=-}f, re = $f? +if’ 
1 2. 2 $. 4J » (40) 
Pee ee ee 
It is seen that rp = 8,1 = & . 


The shearing stress at the wall again differs from the Stokes value only if rs is in- 
cluded; to this approximation, 


r(0) = -(?) i (5) ~ 3 ro) | 41) 


“Yi 1 — 215k"*” (42) 
Tw, 


and we obtain 


whereas for k < 1, from Eq. (37), 


“¥i = 120k. (43) 
Tw, 
As before, both limiting cases are plotted in Fig. 3, together with an estimated smooth 
transition curve. It is seen that the Stokes limit is approached faster by the h-flow than 
by the g-flow; for k = 1 the deviation for the h-flow is about 20%. 
A further set of exact solutions may be obtained if the basic steady, two-dimensional, 
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stagnation-point flow is replaced by a three-dimensional one. It is not necessary to 
assume rotational symmetry; the “potential” flow along the plate may be of the form 
w = az with a, # a, . The steady viscous solution has been obtained by 


“= Q, 2, 
Howarth [5]. If the plate moves or oscillates in any direction in the xz-plane, further exact 
solutions of the Navier-Stokes equation are easily obtained. No examples will be carried 
out for the three-dimensional stagnation point, as the cases treated before are well 
representative of the phenomena that may be expected. 

It is interesting to note that the heat transfer to the plate remains unaffected by 
the motion of the plate in its plane, if the plate temperature is constant. The temperature 
field is obtained as a solution of the equation (in two dimensions) 

oT oT 


+1 


oT _v aT (44) 
ot y Ox 


ig Soe Marla ae 
oy a OY 

where 7 is the temperature and oa is the Prandtl number; dissipation is neglected, as is 
permissible for high temperature differences between the fluid and the plate. For con- 
stant plate temperature (and constant temperature difference between the plate and 
the fluid), the solution with the property 07'/dx = 0 is appropriate. But the plate motion 
under consideration only affects u and leaves v unchanged, so that for 07'/dx = 0 no 
effect on the solution.of Eq. (44) will be felt. 

If unsteady heat transfer is enforced by a variable temperature on the plate, which 
becomes a function of time but not of z, the resultant equation becomes completely 
analogous to that of the h-flow, discussed before. For ¢ = 1, the same solution can be 
used as before; modifications for o ¥ 1 are easily obtained. 

A case which has not been considered so far in this paper is the problem associated 
with a plate oscillating perpendicular to its plane, or the equivalent case of the oscillat- 
ing basic stagnation-point flow. If the f-flow is unsteady, a df’/dt — term appears in 
the non-linear equation (5), so that solutions with a harmonic time-dependence can 
be obtained only to a linearized approximation. Ultimately, however, the unsteady 
part of the flow will influence the development of the steady part of the f-flow, due to 
the non-linearity of Eq. (5). 

The linearized approximation, or the superposition of a small oscillating part to a 
basically steady f-flow has been investigated by Lighthill [6] as a special case of fluctuat- 
ing flow problems with arbitrary velocity distributions. The linearized solution exhibits 
essentially the same features as found in the cases discussed before, inasmuch that a 
quasi-steady type for low frequencies and a high-frequency type approaching Stokes’ 
solution can be distinguished. Lighthill finds as a “limit’’ between these cases, the value 
of k = 5.6 for the reduced frequency, i.e., higher than the limits which have been found 
above for the g- and h-flows. Lighthill also discussed the effects of the time-dependent 
f-flow on heat transfer, which does not vanish, as the v-component of the flow is affected. 

Considering the flow at the stagnation-point of an oscillating airfoil with nose-radius 
R, the value of a is about U/R. For all reduced frequencies of practical interest in 
airfoil flutter, the high frequency approximation is appropriate for all unsteady boundary 
layer phenomena. 

It may be noted that unsteady rigid rotary motion of the plate in its own plane 
leads to problems related to the well-known Kaérman-Cochran case and its generaliza- 
tions; again, however, solutions with harmonic time-dependence can be found only to a 


linearized approximation. 
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STABILITY OF SPHERICAL BUBBLES* 
By GARRETT BIRKHOFF (Harvard University) 


The perturbation equations for a spherical bubble of radius b(¢) are [1, p. 306] 
b(t)bf’ + 3b’()bi — (h — 1)b’’(Ob, = O. (1) 


It is the purpose of this note to give a general stability criterion for the stability of (1). 
Rewriting (1) in the form 


x’ + p(ix’ + q(dz = 0,7 p = 3b’/b, q = —(h — 1)b’’/b, (2 


we consider first the formal identity 


> 


. 4s > "a : 
- {jz + qr} = - ? [2pq + q’], (3) 


which is an easy consequence of (2). 


Tueorem 1. If ¢ < 0, or if g > 0 and 2pq + q’ < 0, then (2) is unstable. If g > 0 
and 2pq + q’ > 0, then (2) is stable. 


Proof. If q(t) < 0, and x(t), x(t) have the same sign for ¢ = ¢ , then they have the 
same sign for all ¢ > ¢, . This is evident since x’(t,) = 0 offers the only possibility for the 
first sign change, and it implies 2’’(t,) = —gz(t,), whence 2’(t, + dt) = — q(t,)x(t,)dt 
has the same sign as x(t, + dt). Hence x(t) grows forever in magnitude; this is of course 
the non-oscillatory case. 

If g(t) > 0, then we are in the oscillatory case. To see this, replace (2) by the self- 


q exp ( [ p at). (4) 


Then we use the Bocher-Priifer variable 0, defined by tan 0 = — Px'/z. Differentiating 
tan 0, using (4), and simplifying, we get 


adjoint 


d(Px’)/dt + Qx = 0, P = exp (/ p at), Q 


.. sin’ 8 > 0. (5) 


d@/dt = Q(t) cos’ 8 + PO: 


*Received Dec. 21, 1954; revised manuscript received March 17, 1955. 
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Hence (if P(¢) is bounded, and Q(t) bounded away from zero), x’ will assume the value 


0 infinitely often, whenever 0 = 7, 27, 37, --- . Now, considering (3), we see that the 
amplitude of successive oscillations increases or decreases as stated in Theorem 1. 

Remarks. Theorem 1 does not seem to be in the literature’, which is mostly con- 
cerned with the special case p = 0. It does not cover cases, like that of the Mathieu 
equation x’ + (a + 6b cos t)x = 0, in which the sign of (2p + q’) oscillates rapidly. In 
the special case that g is constant, then the conditions of Theorem 1 reduce to the familiar 
conditions p > 0 (positive damping) and q > 0 (positive restoring force). 

Finally, in the neighborhood of a regular singular point, where tp(t) = a, + -- 
(g(t) = a. + expanded in MacLaurin series, the condition g > 0 reduces in the 
limit to a, > 0, while 2pq + q’ > O reduces to 2a,(a, — 1)/t* + --- > 0, or (in view of 
a. > 0 andt < 0) toa, < 1. Thus we recapture the stability conditions a, < 1 and 
a, > 0 proved previously [1, p. 308] by use of the indicial equation. In summary, the 
stability conditions of Theorem 1 are a direct generalization of the classical Routh- 


- and 


Hurwitz conditions. 
Now, substituting into Theorem 1 from the case (1)-(2) of a spherical cavity, with 


negligible surface tension and cavity density, we get the stability criteria —b’’/b > 0 and 
(— 6b’b’’ — bb’” + b’b’’)/b? > 0. Since b > 0, these are equivalent to 


b”’ <0, (bb’’’ + 5b’b”) <0; (6) 


the second inequality is equivalent to requiring that (b°b’’) be a decreasing function. 
Evidently, the first inequality is Taylor’s original criterion, that acceleration of the 
interface be towards the lighter fluid. The second inequality expresses the condition 
that there be positive damping’. 
Finally, we consider the general case treated recently by Plesset®. In Plesset’s no- 


tation, b = R and the basic equation involved can be written 


bi’ + (8R’/R)bi — Ab, = 0. (7) 
Substituting in Theorem 1, we get the two stability conditions 
A <0 and 6AR’ + A’R < 0, (8) 


of which the second is equivalent to the condition that R°A be a decreasing function. 
Historical remarks. The stability of spherical bubbles was apparently first con- 
sidered by Riabouchinsky (Proc. Int. Congr. Appl. Mech., Stockholm (1930), p. 149). 
Riabouchinsky suggested that the spherical shape was stable when the bubble pressure 
was less than the ambient pressure, and unstable when it exceeded the ambient pressure. 
Later [Proc. Roy. Soc. A201 192, (1950)], Sir Geoffrey Taylor proposed the more accurate 


1See R. Bellman, Stability theory of differential equations, New York, 1953, and references given there. 
However, C. T. Taam (Proc. Am. Math. Soc. 5, 705-15, Thm. 2 (1954)) has recently obtained a closely 
related result. 

*The fact that negative damping is an important destabilizing factor in cavity collapse was pointed 
out by the author in July, 1952, at the Underwater Ballistics Conference in Pennsylvania State College. 

3M. 8S. Plesset, J. Appl. Phys. 25, 96-8, formula (13) (1954). Plesset’s assertion that his G(t) < 0 
implies stability, seems to be unjustified. A much more elaborate study of bubble instability has recently 
been given by R. H. Pennington, Tech. Res. Rept. 22 (1954), Appl. Math. and Statistics Lab., Stanford 


University. 
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condition that there was stability when b” < 0, and instability when b” > 0. The need 
for a second condition was suggested by the author’ in 1952; in [1], the case of a vapor- 
filled bubble was treated. The present note supplies such a condition for the general 


case that b changes by a large ratio. 
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ERROR BOUNDS FOR A NUMERICAL SOLUTION OF A RECURRING 
LINEAR SYSTEM* 


By A. DE LA GARZA (Carbide and Carbon Chemicals Co., Oak Ridge, Tenn.) 


1. Introduction. Suppose that we are given a bounded region R, f > 0, a function 
¢@, on the boundary B of R, and point functions f and g in R, and that we are required 


to determine ¢ so that 
(P): &o/dx*® + 0°o/dy’ = fo+ginR, @ = o, onB. 


In order to calculate ¢ approximately, we first construct a square grid in a rectangle 
S that has sides parallel to the coordinate axis and is large enough to contain RF in its 
interior. Then we approximate (P) by a system of linear algebraic difference equations, 
arriving at a nonhomogeneous linear algebraic system 


Ag as (1) 
where ¢ is a vector having a component associated with each of the N grid points inside 


R, » is a known N-vector, and A is a known N X N matrix. Finally, we obtain an 
approximate solution ¢» of (1), committing an error 


.= A“'(n — Ady). 


A direct estimate of ¢ is given in Eq. (4) below. Our principal object is to show that 
« can be estimated far more conveniently, though less exactly, from a properly chosen 
system L of linear algebraic equations set up on the M > N grid points interior to the 
rectangle S. Though motivated and illustrated by problem (P), Z can be constructed 
as soon as A (with the properties listed below) is known. This method applies generally 
to any system (1) with those properties. 

2. Estimate for «. We first list properties of A = (a,;): 


a:;<0, t=10N; a, >0, t¥j, i,j = DN; 


Ya;,<0, i= 1N, 


" 
j=1 
the inequality holding for at least one value of 7; 


*Received January 27, 1955. The work described in this paper was done under AEC Contract 
W-7405-eng-26. 
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the matrix A cannot be transformed by the same permutation of rows and columns 
to the form 

PU 

0 Q 


d 
where P and Q are square matrices, and 0 consists of zeros. It follows from Theorem 
II in [1] that | A | 0. 

We now develop a series expression for «. Let m be the diagonal matrix with diagonal 
entries m; = — a;; ,7 = 1(1)N. It may be shown by induction that 


A’ = — >) D’m™' + D"A"," (2) 


where D = I + m™'A. Let A, ,g = 1(1)N, be the characteristic roots of D. From Theorem 
II in [1], it is seen that for | \ | > 1, the determinant | D — J | ¥ 0; therefore, | \, | < 1 
for g = 1(1)N. This implies that D” —- 0 asn — ©. Hence, if ¢o is an approximate solution 


to (1), the error in ¢, being e« = (e,), we have from (1) and the limit in (2) that 
e=- >, D’m~'p, (3) 


where the residual vector p = (r,;) is 7 — Ady. 

We proceed to give a direct estimate for e. Use is made of the abmatrix; see [2]. Let 
B = (b,,) be a matrix. The matrix a(B) = (| b,, |) is called the abmatrix of B. If B = 
(b,,) and C = (c,,), p = 1(1)M, q = 1(1)N, a(B) > a(C) means | b,, | = | ¢,, | for all 
p, q. We see from these definitions that a(A + B) < a(A) + a(B) and a(AB) < a(A) 
a(B), provided the indicated operations are permissible. Finally, when b,, > 0 for all 
p, 9, B = a(B); and the prefix a will not be used. With this notation, since D > 0, 


—_ } 1 . \ 
m; > 0,m; < m_, , and a(p) < Wy rv , Where m, = min(m;), ry = max | r, |, and 
vy = (1,1, --- 1)’ in N dimensions, we have from (3) that 
ale) < (ry/mz) =. D’ by _ (+) 


p=0 


Let us now consider an 1 X M matrix EF > 0, M => N, which contains an N XK N 
submatrix Ez > D and which satisfies E? — 0 as p ©. Let yy = (1,1 ---)’in M 
dimensions. We see that 


= (E’y)* = i D’ by ’ (5) 


where (£’y,;,)* is the column vector formed by the N elements of E”*~y corresponding 
to the N rows of E which contain the submatrix EF, . Since LE? > 0 as p=, (J — EF) 
is non-singular. We may verify that 


> EB’ =(I- EB)". (6) 


p=0 
Hence, 


r= >. Evy (7) 


p=0 
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is the solution of the system 

(I - E)r = Wau ° (8) 
Therefore, if 7* = >°°_, (E’¥y)*, we see from (4), (5), and (6) that 

ale) < r*ry/m,z , (9) 
which is an estimate of the error ¢ in the approximate solution ¢ of (1). In applying 
this bound, we choose / in such a way that (8) has a readily obtainable solution. 

3. Application to a Poisson-type equation. To illustrate use of the methods in part 2, 

we return to the problem (P) in part 1 and its approximation by finite differences. Let 


the grid spacing in S be Ax = Ay = 6. Assign one of the numbers g = 1(1)N as an index 
to each of the grid points interior to R. The first order system approximating (P) then is: 

—m,o, + T,(¢) = g,6 — L,(B), q = 1()N, (10) 
where: 


m, = 4; 


T,(@) is the sum of values of ¢ on non-boundary points adjacent to the gth point; L,(B) 
is a linear form of boundary values required only where the gth point is adjacent to a 
boundary. An example will make this situation clear. With first order linear approxi- 
mation throughout, the equation associated with ¢, in Figure 1 is 


o, — [4 + f.6 + ¢/(1 — ldo + 3 + os = god — o2/(1 — 0). 
ae 


7 
8 
1 








| 


| 

| 
|__| _f\ =a 
tN | | " 


date $3 | | 


| 
| 








/ 











Fia. 1 


We see that the gth row of the matrix D associated with (10) has elements m;' 
corresponding to the unit coefficients of 7',(@), all other elements of D being zero. Let 
the sides of the rectangle S be of length Hé and Jé. Consider the system (11) correspond- 
ing to the (HW — 1) X (J — 1) grid points interior to S: 


bg — me (tecar + bors $+ bce + &-1.0 = 1, i = O(1)/, h = O(1)H, (11) 
tio = thor = to .¢ = ts = 0, 


where m, = min(m,). Write the system (11) in the matrix form (8), the equations 
corresponding to the grid points interior to R being written in the same order as in (10). 
We see that the resulting matrix E then contains an N X N submatrix Ep, > 0 with 
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entries mz,‘ in the same positions as the entries m,' of the matrix D associated with 
(10). Hence, Ez > D. From Theorem II in [1], we see that the characteristic roots of E 
are less than one in absolute value; hence, E’ — 0 as p —~-. It follows that the solution 
of (11) may be used in (9) to compute an error bound for an approximate solution of 
(1). We may verify that the solution of (11) is 

r=H-1,s=J-1 


t. =(HD" > — C,,sin hrv/H sin isr/I, (12) 


er Eo 
where 

C.,, = {1 -— (-)D’'][1 — (-—)D‘]smrr/H sin sr/T} 

-{(1 — cosrm/H)(1 — cos sx/I)[1 — 2mz'(cosra/H + cos sr/I)]}"". 


For specifying the largest absolute residual which may be tolerated in an iterated 


solution, a bound for the maximum error ey = max | e; | is useful. From (9), 
er < ty? mr ; (13) 
where f, = max(t¢,.;). For simplicity, let both H and J be even. We may verify that 


t, occurs ath = H/2,i = I/2, and that 


r=H-1,s=I] 


ty = —4(HI)™ 2. (—1)°*?11 — 2mz'(cos rr/H + cos sx/I)}"' 


-cot rr/2H cot sx/2/, (14) 


where the summation is on odd values of r and s. For illustration, ty in (14) with m, = 4 
is tabulated in Table 1 for H, J = 8 (2) 20. Since ¢y in (14) increases as m,, decreases, 
the tabulated values may be used in (13) to state error bounds for an approximate 
solution of the linear system (10). 


TABLE 1 
Factor ty for use in (13) 
H, I = 8(2) 20 


H/I 8 10 12 14 16 18 20 
8 18.6 

10 23.7 29.2 

12 25.6 34.5 42.2 

14 27.6 38.5 48.6 B75 

16 29.0 41.5 53.8 65.1 75.2 

18 30.0 43.8 57.9 71.4 84.0 95.2 

20 30.6 45.4 61.1 76.6 91.4 105 118 


Acknowledgements. Mr. E. B. Carter of the Numerical Analysis Department at 
K-25 computed Table 1. The referee’s comments simplified and generalized the original 


results. 
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SOUND TRANSMISSION FROM A TUBE WITH FLOW* 
By G. F. CARRIER (Harvard University) 


1. Introduction. The increasing technological importance of thermo-acoustic 
phenomena in combustion chambers and other apparatus leads naturally to the need 
for understanding the transmission and reflection of acoustic waves at the inlet and 
exit sections of tubes through which a gas is flowing at moderate Mach number. In 
this paper, we treat this question for the inviscid perfect gas. The flow fields adopted 
are the simplest which are reasonable approximations to interesting physical situations 
and which lead to tractable problems. The problems are resolved by an extension of 
the work of Schwinger and Levine [1] (their work is essentially the M = 0 problem), 
and, in some interesting cases, our results can be expressed directly in terms of theirs. 
The modifications in these cases are reminiscent of subsonic aerodynamic results ob- 
tained using the Prandtl-Glauert method. 

2. Formulation of the problems. The configuration to be treated here is depicted 
in Fig. (2.1). The steady flow through and about the tube is characterized by the specifi- 


" 
r 














Fic. 2.1. Geometry of sound transmission problem. 


cation of the thermodynamic state variables and the velocity field; the velocity field, 
in turn, is characterized by a single non-vanishing velocity component in the z, direction. 
This velocity has the value v, for r, < R and v, for r; > R. In two of the problems of 
interest v, = v, and the state variables are the same in each region. In the third v, > v, 
and the state variables may differ in the two domains. The first two cases are distinguished 
according to whether v, is greater or less than zero. 

If we restrict our analysis to small disturbances, assuming a perfect, inviscid, non- 
heat conducting gas, the equations governing the propagation of the disturbance in 
each region take the well known form 


p div v’ + pi + vp!, = 0, (2.1) 
(2.2 


I 
= 


p(vi + vv’) + grad p’ 


*Received March 11, 1955. 
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Here, the primed quantities (v’, p’, p’) are the perturbations in velocity, density and 
pressure, p is the steady state (constant) density, and the isentropic pressure-density 
law implies that p’ = a’p’ where a is the acoustic speed associated with the steady 
flow state; z, , 7: are the space coordinates. Each of these quantities will require the 
subscript 1 or 2 to distinguish the domain to which it applies when we discuss the flow 
in which v, ¥ 2. 

It is convenient to introduce the following notation: v’ = aR’ grad ¢, R’ = 
Ri — M’)'”,s = 2,/R’,y = 7,/R, k’ = w R*/a'(1 — M”), and¢ = Y(s, y) exp { —tw[t + 


Mz,/a(1 — M”’)]}. It follows from these and the foregoing that 
Ay+ky=0 (2.3) 
and 
p’ = —pa (My, — iky) exp {—iw[t + Mz,/a(1 — M”)]}, (2.4) 


where A is the cylindrical coordinate Laplace operator." The boundary conditions we 
shall use differ in the three problems. If problem I is characterized by v, = v, > 0, 
the essential information is: p’(z, R, t) is continuous on s > 0, y,(s, 1) = Oons < 0. 
Note that, as in oscillating airfoil theory, one can anticipate that y(s, 1) will not be 
continuous for s > 0. In problem II, v, = v, < 0 and both p’ and y are continuous on 
y = 1, s > 0; again, y,(s, 1) = 0 whens < 0. For problem III, with v, and v, positive 
but not equal, the boundary conditions are the same as those of problem I. 

3. The continuous flow problems. Problems I and II differ only in minor respects 
from the problem treated in [1]. This enables us to reduce these problems to such a 
form that the results of interest can be obtained directly in terms of quantities computed 
in that paper. The ordinary differential equation which x, the Fourier transform of y, 


lie. x(—, y) = f2. ¥(a, y) exp (—ztx) dx], must obey is 
y (yx), + (kh -— &)x = 0. (3.1) 
The appropriate solution (obeying the radiation condition) is 
x = Alé)Jo([k? — &]'’’y), y¥<l 
and 
x = B)Hy ((k’ — #)'y), o> i. (3.2) 
The boundary conditions can be written concisely if we define h(é) = (tk — «M#) 
Ix(é, 1—) — x(é, 1+)]. The pressure boundary condition then merely implies that 


h(&) is the transform of a function, H(x), which is zero for x > 0. We also note that 
w(t) = x,(é, 1) is the transform of a function, W(x) which vanishes for x < 0. Now 


2 


h(t) = i(k — Mt){A®J([F — £]'”) — BOHS ([k — £]'”)} 


and 
w(t) = —A(h — 870, ((kP — &)'”) = -—BR — &)'7H) (kh -— &]'”). 


Eliminating A and B, we obtain (using J,(¢)Ho(¢) — Jo(OHi? (6) = 2i/ro), 
—ri(k? — &)A((k — &)) J.((k — &)' he) = 2i(k — Mé)w(€). (3.3) 
Our foregoing remarks imply that h(é) is analytic when Im £ > 0 and that w(£) is analytic 


1We consider only problems in which the incident wave is a plane wave in the tube normal to its 


axis and approaching the origin. 
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when Im ~ < 0. The domains of analyticity of these functions can be extended if we 
associate with k a complex value (Im k = e > 0). When this is done, we can anticipate 
that h, w, and the coefficient of h in (3.3) will be analytic in the strip | Im ¢ | < e. This 
is demonstrated in detail in [1]. We now write’ 

L(g) = wily” ((k? — #)') Jit’ — &]'”) = L.@/L-©, 
where L, and L_ are respectively analytic above —ie and below ze. One can now apply 


the conventional arguments of the Wiener-Hopf technique to find h(é), and these argu- 
ments lead to the conclusion that each side of the equation 


~(k? — E)L.(OhE) = 2i(k — ML_)w() (3.4) 


is the analytic continuation of the other, and that each represents an entire function of &. 
When W > 0 (the exit problem), it is required that H(x) — 0 as x — 0— and we 
can anticipate that H(x) ~ | 2 |'” in this neighborhood. It was shown in [1] that 
L.(g) ~ &"” for large — (with Im £ > 0). It follows that the entire function, ( — 
k°)L.(£)h(£), is a constant (because it is bounded at ~). 
When MV < 0, the situation is different. In this case, instead of demanding that 


h(x) — 0 as x — O—, we require that y be continuous on the extension of the tube. 
This implies that —7zh/(k — Mé) be the transform of a function which behaves like 
z'’* near « = 0 and, in this case, the entire function (k° — £)L,()h(é), is a constant 
multiple of (k — Mé). One also can see this by noting that, in the entrance case, the 

1/2 


radial velocity perturbation near the edge must behave like 2~'’*; then, since L_(¢) ~ & 
(see [1]), the foregoing result ensues. 

We are primarily interested in the reflection coefficient associated with the tube end. 
Since the pressure field in the tube for large negative x must be of the form H(x) ~ 
exp (ks) + N exp (—iks), it follows that the reflection coefficient N is the ratio of the 
residues of h(£) exp (7ts) at —k and +k, respectively (ie. N = —L,(k)/L.(—k) for 
the exit problem and N = —[(1 + )/(1 — M)JL,(k)/L.(—k) for the inlet problem, 
where ./ < 0). For the exit problem, these residues are precisely those computed in [1]. 
Thus our V{(wR/a)(1 — M)’”] is precisely the R of [1]. Note that this is the reflection 
coefficient for the pressure field (and, as it happens, for the velocity field, too) but not 
for the potential. 

For the inlet flow, on the other hand, the residues of h(€) exp (és) at +k are not 
the same as those of the exit flow but rather those of this expression multiplied by 
i(k — Me). Thus, for WZ < 0, the reflection coefficient N’(k) is given by [(1 + M)/(1 — 
M)\N(i 

It is interesting to note that the reflection at the inlet end of the tube is much poorer 
than that at the exit end. This is most readily rationalized by noting that, at the inlet 
end, the incident wave “sees” a much smaller wave length to radius ratio. This is not 
an adequate basis for an estimate, however, because the non-convective theory gives 
a quadratic dependence on this ratio for small frequency whereas here the dependence 
is stronger and, in particular, N is not unity for w = 0. 

The major result is to notice that for low frequencies, a moving tube has much 
greater radiation losses from an inlet end than it does in a stationary situation but that 
the losses at the exit are only affected when M gets reasonably close to unity. 


2The notation is nearly identical with that of [1] in order that the reader may readily translate the 
thorough discussion of that paper into the straightforward generalization of this section. 
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Figures (3.1) and (3.2) are graphs of N and 1/R where N = —| N | exp (27kl/R). 
These are taken from [1]. 
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4. The jet flow. In this section we consider the algebraically more difficult problem 
where M, + M, and a, + a, . In this problem, it is more convenient to distort the coordi- 
nate system uniformly (i.e. we use the same scale for x and r) and we introduce the nota- 
tion: s = 2/R, y = r/R, m; = 0;/a; , k; = wR/a; , Bj = 1 — Mj andv; = a; grad 
¢;(s, y) exp (—iwt). The differential equations governing the ¢; become (proceeding 
from Eqs. (2.1) and (2.2)) 


Ad, — (a1, = - ity) ¢, = 0, (4.1) 
os 


5 _\2 | 
Ad, — (1, z= _ ity) gd. = 0. (4.2) 


2 again introduce the Fourier transforms with regard to x o e @¢;, ca 1e 
If we again introduce the F t f th regard t f the ¢; ll them 
x; (&, y), and integrate the ordinary differential equations which govern the x; (&, y), 
we obtain 


xilé, y) = Alé)Jo(e.y), (4.3) 
xo(é, y) = B(é)Hy’ (zy), (4.4) 

where z; = (k; — M,&)? — #. The velocity at y = 1, W(s) has the transform 
wt) = —z,AJ,(z;) = —z,BH{"’ (22) (4.5) 


and the jump’ in pressure at y = 1, H(s), has the transform 
h(t) = (k, — M,&)AJ,(z:) — (kz — Mt)BH%" (z,). (4.6) 
Equations (4.5) and (4.6) may be combined to eliminate A and B and to obtain 
ziL(é)h(—) = w(é)[a — (M, + M.£,/8.)€], (4.7) 


*This is the dimensionless jump (p{ — p$)/pi. Note that p: = p2 for equilibrium of the steady flow. 
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where 

fae — (My + M8, /8,)ENo/20) Fle) (4.8) 
(k, — M,£)z,H§? (22) J (2:) — (ki — Mié)2z2H}" (22) J o(2:) 

and where a is so chosen that L/(é) is finite at that zero of the denominator of Eq. (4.8) 


which, when states 1 and 2 are identical, is located at & = k./M, . 
It is convenient to factor L(é) in the following manner. 





L(é) 


L(t) = wid ,(z,;)H{"(z,), L*(&) = L/L, , 


—t[a ~ (M, + M,£, ‘B2) Jeo" (22) — ——— (4 9) 
r[(k, — M,8)z,Jo(z,)H{" (e:)2oH\” (22) — (ke — Mot)H5" (2.)21H}” (2:)J(2:)] ? 


and the factor L,(&) is precisely the Schwinger-Levine L(¢) where £ = 6, + hy . 

The solution of Eq. (4.7) for h(é) now requires that we split L(¢) [and hence both 
L,(é) and L*(£)] into factors L,(&), L_(é) which are analytic above —ie, and below ie, , 
respectively (here e; is the imaginary part of the root of z; = 0 when Im k, = e). 

The reflection coefficient will be given again by the ratio of the amplitudes of the 
returning and the outgoing waves in the tube. The inversion integral for the pressure 
H(x) in the tube (as x — — ©, the external pressure tends to 0 so the interior pressure 


becomes the pressure jump H(2)) is given by 


l . 
H(2) = 5 [ (L.@Y" exp (a) le@? ae 
and the asymptotic result for —2 >> 1 is given by the residues at the zeros of z} . That is, 
the reflection coefficient is again L,(—b.)/L.(b,) where —b, and b, are the zeros of 2; . 
However, L. is composed of two factors, the L.(£) associated with ZL, as used in the 
foregoing and the factor L* . The Cauchy integral formula gives L* in the form 


In LX) = = [ (€ — )"' n L*@®) a (4.10) 
aT oo 


and, in particular, we define the quantity Q so that 


In Q = In(L*(—b,)/L*(b,)) = —t [ 2 (4.11) 

. miB; Jw (—E — bE + b.) ~”’ 

where the path of integration is indented to pass below —b, and b, . When k, = k, , 

M, = M, , L*(é) is identically unity and the contribution of (4.11) is Q = 1. For other 

cases, however, L*(£) — 1 as  ~ +o and the convergence of the integral is assured. 
The reflection coefficient in such cases is given by the product 


Qk: , ke , My , M.)N(k,/6,). 


Again one can use the results of [1] to find N but the determination of Q requires the 
numerical evaluation of the integral of (4.11) for those values of the four independent 
parameters a, , a, , M, , M, which are of interest. 
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SOLUTIONS OF THE HYPER-BESSEL EQUATION* 
By CHIA-SHUN YIH, (State University of Iowa) 


In problems of hydrodynamic stability involving axial symmetry, it is sometimes 
necessary to find the solutions of a differential equation of the type 


in which n is a positive integer, and (with D = d/dr) 
L=D?+r'D-r?*-» 
is the Bessel operator of the first order. In this note, solutions of the equation 
Lsf = 0 (1) 


in which 


L=D+r'D-pr’+Fk 


will be given explicitly. The theorem one seeks to establish is the following: If p (taken 
to be positive for convenience) is not an integer, the solutions of Eq. (1) are r"J (54m) (Kr) 
in which m = 0, 1, 2, --- , m — 1; otherwise they are r”J,,,,(kr) and r"N,,.,,(kr), with 
m ranging over the same integers. The symbols J and N stand for the Bessel function 
and the Neumann function, respectively. 

Proof: It is known that the solutions of L,f = 0 are J.,(kr) for p not equal to an 
integer and J,(kr) and N,(kr) for p equal to an integer. Thus it suffices to show that 
if r°Z,..(kr) (in which Z stands for either J or N) satisfies L;*'f = 0, then r°“'Z,.,.,(kr) 
satisfies L3*’f = 0, since the proof for r”J_;,m) (kr) is identical with that for r”Z,.,,(kr). 
This will be accomplished if one can show that L,r**'Z,,.,.,(kr) is equal to a constant 
times 7°Z,.,(kr). By straightforward differentiation one has 


Lor Zosesikr) = 2" LD Zoysacilkr) + 8(8 + 1)r’'Z,4041(kr) 
+ 2(s + 1)r’°DZ,.04:(kr) + (8 + Ir? "Z,404i(kr) 
= 9D Zoseailkr) + (8 + 1)'r°"Z,404:(kr) 
+ 2(s + 1)r’DZ,.,.,:(kr). 
But 


and [1] 


DZ o.1+\(kr) = x ESet lg aay 4 Zot) |. 


kr 





*Received May 14, 1955. 
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So 


Ly Zosesikr) = 1°" "Lyser:Zoresi(kr) + [2p(s + 1) + + 1) 'Z, 4 44:(kr) 


“ 


+ (s + 1)r"'Z,404:(kr) + 2s + Dri - poet Zorerilkr) + Zk) 


= 2(s + 1)kr*Z,,.(kr) 
since 
| SOE | = 0 
by definition of Z. 
Dr. Y. C. Fung of the California Institute of Technology communicated to the 
writer a different proof of the present result by means of Almansi’s theorem [2] on hyper- 


harmonic functions. His proof will not be presented here. 
It may be noted that since [1] 


9 
Z,-(kr) + Zy.i(kr) = a Z,(kr) (2) 

and since by the theorem just proved rZ,,,(kr) and Z,(kr) are solutions of 
Lif = 0, (3) 


it follows from Eq. (2) that rZ,_,(kr) is also a solution of Eq. (3). In fact, by repeated 
use of Eq. (2) and a similar one obtained by changing p to — pin Eq. (2), it can be proved 
that if the m in the subscripts of the solutions given in the theorem is changed to —m, 
the results will still be solutions of Eq. (1). These solutions are of course not independent 
of the ones given in the statement of the theorem. 
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BOOK REVIEWS 
(Continued from p. 392) 


Mathematics of engineering systems. By Derek F. Lawden. Methuen & Co., Ltd., London, 
and John Wiley & Sons, Ine., New York, 1954. viii + 380 pp. $5.75. 


This book provides a course in applied mathematics in which a selection of topics of interest in 
modern engineering and applied science are treated on an intermediate level. After an introductory 
chapter which includes complex numbers there are two chapters on linear ordinary differential equations 
with constant coefficients. The classical methods are treated in the first of these, with applications to 
problems of servomechanisms and stability; modern methods are taken up in the second, involving the 
use of standard inputs (unit steps, unit impulse, and sinusoidal functions) and the Laplace transform. 
Chapter 4 is concerned with Fourier series and integrals and applications of Fourier transforms. The 
final chapter provides a brief introduction to non-linear differential equations, in which main attention 
is paid to Van der Pol’s and to Duffing’s equation. This is a book of mathematics, not of engineering; 
but the large number of examples refer to practical engineering problems. They are drawn mostly from 
electrical engineering (electronic amplifiers and oscillators, electric circuits, and servomechanisms), but 


the book should be of interest to engineers and scientists in many fields. 
P. S. Symonps 


La Théorie des Fonctions de Bessel. By Gerard Pétiau. Renseignements et Vente au 
55. 477 pp. $7.20. 


Service des Publications du CNRS, Paris, 1! 


This book would be a very valuable addition to the literature on Bessel functions if it were not 
lacking some rather vital features. The most serious omission is an alphabetical index. This is seldom of 
much value in trying to find an equation but it does serve an important function in finding some things 
and the lack of one seriously reduces the usefulness of the book. The book contains references for most 
topics of recent development and for certain topics of related interest to the main topics but there is no 
general list of references and by and large the system of references is rather weak. The book is also 
lacking a preface; thus there is nowhere in the book any discussion of the scope or extent to which this 
book supersedes previous works on Bessel functions or even for what purpose the book was intended. 

On the more favorable side, the book does contain a tremendous collection of formulas. Most of 
the derivations are short and the book seems to be arranged with this particularly in mind. Many formulas 
and even some interesting topics are not to be found in some of the standard references (for example 
Watson’s “The Theory of Bessel Functions” or the recent Bateman Manuscript Project ‘Higher 
Transcendental Functions’) There is a short chapter on Bessel integral functions 

Ji,(z) = — | u'J,(u) du 
and related integrals. There are also some definite integrals that are not easily found elsewhere. 

The scope of the book is extensive. The first 328 pages are devoted to the properties of Bessel functions 
and related functions. It includes much of the material in Watson’s treatise plus some more recent topics. 
The next 107 pages describes some of the principal applications particularly in relation to the solution 
of partial differential equations. By comparison this part is much more illustrative than exhaustive. 
In addition there are 14 pages of short tables and 11 pages of graphs. 

The derivations of some of the formulas employ standard mathematical techniques which are of a 
somewhat advanced nature. This, coupled with the scarcity of references, will make the book difficult 


to follow for a person not well trained in mathematics. 
G. F. NEWELL 


























